1module costFunction_mod
2 ! MODULE costFunction_mod, (prefix='cfn' category='5. Observation operators')
3 !
4 !:Purpose: To compute Jo term.
5 !
6 use midasMpi_mod
7 use obsSpaceData_mod
8 use rttov_const, only : inst_name, platform_name
9 use tovsNL_mod
10 use utilities_mod
11 use obserrors_mod
12 use codtyp_mod
13
14 implicit none
15
16 save
17 private
18
19 ! public procedures
20
21 public :: cfn_calcJo, cfn_sumJo
22
23 integer, allocatable :: channelNumberList(:,:)
24 character(len=15), allocatable :: sensorNameList(:)
25
26contains
27
28 !--------------------------------------------------------------------------
29 ! cfn_calcJo
30 !--------------------------------------------------------------------------
31 subroutine cfn_calcJo(lobsSpaceData)
32 !
33 !:Purpose: To compute JO contribution of each assimilated and diagnosed
34 ! datum, and to store the result in OBS_JOBS
35 implicit none
36
37 ! Arguments:
38 type(struct_obs), intent(inout) :: lobsSpaceData
39
40 ! Locals:
41 integer :: bodyIndex
42
43 !$OMP PARALLEL DO PRIVATE(bodyIndex)
44 do bodyIndex=1,obs_numbody(lobsSpaceData)
45
46 if ( obs_bodyElem_i( lobsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated) then
47 call obs_bodySet_r( lobsSpaceData, OBS_JOBS, bodyIndex, &
48 ( obs_bodyElem_r( lobsSpaceData, OBS_WORK, bodyIndex ) &
49 * obs_bodyElem_r( lobsSpaceData, OBS_WORK, bodyIndex ) &
50 ) / 2.d0 )
51 else
52 call obs_bodySet_r(lobsSpaceData,OBS_JOBS,bodyIndex, 0.0d0)
53 end if
54
55 end do
56 !$OMP END PARALLEL DO
57
58 end subroutine cfn_calcJo
59
60 !--------------------------------------------------------------------------
61 ! cfn_sumJo
62 !--------------------------------------------------------------------------
63 subroutine cfn_sumJo( lobsSpaceData, pjo, beSilent_opt )
64 !
65 !:Purpose: To compute the sum of Jo contributions saved in OBS_JOBS. Also,
66 ! to compute contribution of each family of observation (for
67 ! diagnostic purposes)
68 implicit none
69
70 ! Arguments:
71 type(struct_obs), intent(in) :: lobsSpaceData
72 real(8), intent(out) :: pjo ! Total observation cost function
73 logical, optional, intent(in) :: beSilent_opt
74
75 ! Locals:
76 integer :: bodyIndex, tovsIndex, sensorIndex, headerIndex, bodyIndexBeg, bodyIndexEnd
77 integer :: channelNumber, channelNumberIndexInListFound, channelIndex
78 integer :: sensorIndexInList, sensorIndexInListFound
79 logical :: beSilent
80
81 real(8) :: dljoraob, dljoairep, dljosatwind, dljoscat, dljosurfc, dljotov, dljosst, dljoice
82 real(8) :: dljoprof, dljogpsro, dljogpsztd, dljochm, pjo_1, dljoaladin, dljohydro, dljoradar
83 real(8) :: dljotov_sensors( tvs_nsensors )
84 real(8) :: joTovsPerChannelSensor(tvs_maxNumberOfChannels,tvs_nsensors)
85 character(len=15) :: lowerCaseName
86
87 logical :: printJoTovsPerChannelSensor
88
89 real(8), allocatable :: joSSTInstrument(:)
90 integer, allocatable :: nobsInstrument(:), nobsInstrumentGlob(:)
91 integer :: SSTdatasetIndex, codeType, ierr
92
93 if ( present(beSilent_opt) ) then
94 beSilent = beSilent_opt
95 else
96 beSilent = .false.
97 end if
98
99 if (.not. allocated(channelNumberList)) then
100 allocate(channelNumberList(tvs_maxNumberOfChannels,tvs_nsensors))
101 end if
102 if (.not. allocated(sensorNameList)) then
103 allocate(sensorNameList(tvs_nsensors))
104 end if
105
106 if(oer_getSSTdataParam_int('numberSSTDatasets') > 0) then
107 allocate(joSSTInstrument(oer_getSSTdataParam_int('numberSSTDatasets')))
108 allocate(nobsInstrument(oer_getSSTdataParam_int('numberSSTDatasets')))
109 allocate(nobsInstrumentGlob(oer_getSSTdataParam_int('numberSSTDatasets')))
110 end if
111
112 call readNameList
113 printJoTovsPerChannelSensor = .false.
114 if (any(sensorNameList(:) /= '') .and. any(channelNumberList(:,:) /= 0)) then
115 printJoTovsPerChannelSensor = .true.
116 end if
117
118 dljogpsztd = 0.d0
119 dljoraob = 0.d0
120 dljoairep = 0.d0
121 dljosatwind = 0.d0
122 dljosurfc = 0.d0
123 dljoscat = 0.d0
124 dljotov = 0.d0
125 dljogpsro = 0.d0
126 dljoprof = 0.d0
127 dljochm = 0.d0
128 dljosst = 0.0d0
129 dljoaladin = 0.d0
130 dljoice = 0.0d0
131 dljotov_sensors(:) = 0.d0
132 joTovsPerChannelSensor(:,:) = 0.0d0
133 dljohydro = 0.0d0
134 dljoradar = 0.0d0
135 joSSTInstrument(:) = 0.0d0
136 nobsInstrumentGlob(:) = 0
137 nobsInstrument(:) = 0
138
139 pjo = 0.0d0
140
141 do bodyIndex = 1, obs_numbody(lobsSpaceData)
142
143 pjo_1 = obs_bodyElem_r(lobsSpaceData, OBS_JOBS, bodyIndex)
144
145 ! total observation cost function
146 pjo = pjo + pjo_1
147
148 ! subcomponents of observation cost function (diagnostic only)
149 select case(obs_getFamily(lobsSpaceData, bodyIndex_opt = bodyIndex))
150 case('UA')
151 dljoraob = dljoraob + pjo_1
152 case('AI')
153 dljoairep = dljoairep + pjo_1
154 case('SW')
155 dljosatwind = dljosatwind + pjo_1
156 case('SF')
157 dljosurfc = dljosurfc + pjo_1
158 case('SC')
159 dljoscat = dljoscat + pjo_1
160 case('TO')
161 dljotov = dljotov + pjo_1
162 case('RO')
163 dljogpsro = dljogpsro + pjo_1
164 case('PR')
165 dljoprof = dljoprof + pjo_1
166 case('GP')
167 dljogpsztd = dljogpsztd + pjo_1
168 case('CH')
169 dljochm = dljochm + pjo_1
170 case('TM')
171 dljosst = dljosst + pjo_1
172 case('AL')
173 dljoaladin = dljoaladin + pjo_1
174 case('GL')
175 dljoice = dljoice + pjo_1
176 case('HY')
177 dljohydro = dljohydro + pjo_1
178 case('RA')
179 dljoradar = dljoradar + pjo_1
180 end select
181 enddo
182
183 do tovsIndex = 1, tvs_nobtov
184 headerIndex = tvs_headerIndex( tovsIndex )
185 if ( headerIndex > 0 ) then
186 bodyIndexBeg = obs_headElem_i(lobsSpaceData, OBS_RLN, headerIndex)
187 bodyIndexEnd = obs_headElem_i(lobsSpaceData, OBS_NLV, headerIndex) + bodyIndexBeg - 1
188 sensorIndex = tvs_lsensor (tovsIndex)
189
190 sensorIndexInListFound = 0
191 if ( printJoTovsPerChannelSensor ) then
192 loopSensor1: do sensorIndexInList = 1, tvs_nsensors
193 call up2low(sensorNameList(sensorIndexInList),lowerCaseName)
194
195 if ( trim(lowerCaseName) == trim(inst_name(tvs_instruments(sensorIndex))) ) then
196 sensorIndexInListFound = sensorIndexInList
197 exit loopSensor1
198 end if
199
200 end do loopSensor1
201 end if
202
203 do bodyIndex = bodyIndexBeg, bodyIndexEnd
204 pjo_1 = obs_bodyElem_r(lobsSpaceData, OBS_JOBS, bodyIndex)
205 dljotov_sensors(sensorIndex) = dljotov_sensors(sensorIndex) + pjo_1
206
207 if ( printJoTovsPerChannelSensor .and. &
208 sensorIndexInListFound > 0 ) then
209 call tvs_getChannelNumIndexFromPPP( lobsSpaceData, headerIndex, bodyIndex, &
210 channelNumber, channelIndex )
211 channelNumberIndexInListFound = utl_findloc(channelNumberList(:,sensorIndexInListFound), &
212 channelNumber)
213
214 if ( channelNumberIndexInListFound > 0 ) then
215 joTovsPerChannelSensor(channelNumberIndexInListFound,sensorIndexInListFound) = &
216 joTovsPerChannelSensor(channelNumberIndexInListFound,sensorIndexInListFound) + &
217 pjo_1
218 end if
219
220 end if
221
222 end do
223 end if
224 end do
225
226 if(oer_getSSTdataParam_int('numberSSTDatasets') > 0) then
227 do headerIndex = 1, obs_numheader(lobsSpaceData)
228 codeType = obs_headElem_i(lobsSpaceData, OBS_ITY, headerIndex)
229 bodyIndexBeg = obs_headElem_i(lobsSpaceData, OBS_RLN, headerIndex)
230 bodyIndexEnd = obs_headElem_i(lobsSpaceData, OBS_NLV, headerIndex) + bodyIndexBeg - 1
231 do bodyIndex = bodyIndexBeg, bodyIndexEnd
232 pjo_1 = obs_bodyElem_r(lobsSpaceData, OBS_JOBS, bodyIndex)
233 dataset_loop: do SSTdatasetIndex = 1, oer_getSSTdataParam_int('numberSSTDatasets')
234 if (codeType == oer_getSSTdataParam_int('codeType', SSTdatasetIndex) .and. codeType /= codtyp_get_codtyp('satob')) then
235 joSSTInstrument(SSTdatasetIndex) = joSSTInstrument(SSTdatasetIndex) + pjo_1
236 nobsInstrument(SSTdatasetIndex) = nobsInstrument(SSTdatasetIndex) + 1
237 exit dataset_loop
238 else
239 if (obs_elem_c(lobsSpaceData, 'STID', headerIndex) == oer_getSSTdataParam_char('sensor', SSTdatasetIndex)) then
240 joSSTInstrument(SSTdatasetIndex) = joSSTInstrument(SSTdatasetIndex) + pjo_1
241 nobsInstrument(SSTdatasetIndex) = nobsInstrument(SSTdatasetIndex) + 1
242 exit dataset_loop
243 end if
244 end if
245 end do dataset_loop
246 end do
247 end do
248 end if
249
250 call mmpi_allreduce_sumreal8scalar( pjo, "GRID" )
251 call mmpi_allreduce_sumreal8scalar( dljoraob, "GRID" )
252 call mmpi_allreduce_sumreal8scalar( dljoairep, "GRID" )
253 call mmpi_allreduce_sumreal8scalar( dljosatwind, "GRID" )
254 call mmpi_allreduce_sumreal8scalar( dljosurfc, "GRID" )
255 call mmpi_allreduce_sumreal8scalar( dljoscat, "GRID" )
256 call mmpi_allreduce_sumreal8scalar( dljotov, "GRID" )
257 call mmpi_allreduce_sumreal8scalar( dljogpsro, "GRID" )
258 call mmpi_allreduce_sumreal8scalar( dljoprof, "GRID" )
259 call mmpi_allreduce_sumreal8scalar( dljogpsztd, "GRID" )
260 call mmpi_allreduce_sumreal8scalar( dljochm, "GRID" )
261 call mmpi_allreduce_sumreal8scalar( dljosst, "GRID" )
262 call mmpi_allreduce_sumreal8scalar( dljoaladin,"GRID")
263 call mmpi_allreduce_sumreal8scalar( dljoice, "GRID" )
264 call mmpi_allreduce_sumreal8scalar( dljohydro, "GRID" )
265 call mmpi_allreduce_sumreal8scalar( dljoradar, "GRID" )
266 do sensorIndex = 1, tvs_nsensors
267 call mmpi_allreduce_sumreal8scalar(dljotov_sensors(sensorIndex), "GRID")
268 end do
269 if (printJoTovsPerChannelSensor) then
270 loopSensor2: do sensorIndex = 1, tvs_nsensors
271 if (trim(sensorNameList(sensorIndex)) == '') cycle loopSensor2
272
273 call mmpi_allreduce_sumR8_1d(joTovsPerChannelSensor(:,sensorIndex), "GRID")
274 end do loopSensor2
275 end if
276
277 ! SST data per instrument
278 do SSTdatasetIndex = 1, oer_getSSTdataParam_int('numberSSTDatasets')
279 call mmpi_allreduce_sumreal8scalar(joSSTInstrument(SSTdatasetIndex), "grid")
280 call rpn_comm_allreduce(nobsInstrument(SSTdatasetIndex), nobsInstrumentGlob(SSTdatasetIndex), &
281 1, "mpi_integer", "mpi_sum", "grid", ierr)
282 end do
283
284 if ( mmpi_myid == 0 .and. .not. beSilent ) then
285 write(*,'(a15,f30.17)') 'Jo(UA) = ', dljoraob
286 write(*,'(a15,f30.17)') 'Jo(AI) = ', dljoairep
287 write(*,'(a15,f30.17)') 'Jo(SF) = ', dljosurfc
288 write(*,'(a15,f30.17)') 'Jo(SC) = ', dljoscat
289 write(*,'(a15,f30.17)') 'Jo(TO) = ', dljotov
290 write(*,'(a15,f30.17)') 'Jo(SW) = ', dljosatwind
291 write(*,'(a15,f30.17)') 'Jo(PR) = ', dljoprof
292 write(*,'(a15,f30.17)') 'Jo(RO) = ', dljogpsro
293 write(*,'(a15,f30.17)') 'Jo(GP) = ', dljogpsztd
294 write(*,'(a15,f30.17)') 'Jo(CH) = ', dljochm
295 write(*,'(a15,f30.17)') 'Jo(TM) = ', dljosst
296 write(*,'(a15,f30.17)') 'Jo(AL) = ', dljoaladin
297 write(*,'(a15,f30.17)') 'Jo(GL) = ', dljoice
298 write(*,'(a15,f30.17)') 'Jo(HY) = ', dljohydro
299 write(*,'(a15,f30.17)') 'Jo(RA) = ', dljoradar
300 write(*,*) ' '
301 if ( tvs_nsensors > 0 ) then
302 write(*,'(1x,a)') 'For TOVS decomposition by sensor:'
303 write(*,'(1x,a)') '# plt sat ins Jo'
304 do sensorIndex = 1, tvs_nsensors
305 write(*,'(i2,1x,a,1x,a,1x,i2,1x,f30.17)') sensorIndex, inst_name(tvs_instruments(sensorIndex)), &
306 platform_name(tvs_platforms(sensorIndex)), &
307 tvs_satellites(sensorIndex), &
308 dljotov_sensors(sensorIndex)
309 end do
310 write(*,*) ' '
311 end if
312
313 ! print per channel information
314 if ( tvs_nsensors > 0 .and. printJoTovsPerChannelSensor ) then
315 write(*,'(1x,a)') 'For TOVS decomposition by sensor/channel:'
316 write(*,'(1x,a)') 'index sensorName channel Jo'
317 loopSensor3: do sensorIndex = 1, tvs_nsensors
318 if ( trim(sensorNameList(sensorIndex)) == '' ) cycle loopSensor3
319
320 loopChannel: do channelIndex = 1, tvs_maxNumberOfChannels
321 if ( channelNumberList(channelIndex,sensorIndex) == 0 ) cycle loopChannel
322
323 write(*,'(i2,1x,a,1x,i4,1x,f30.17)') sensorIndex, &
324 sensorNameList(sensorIndex), &
325 channelNumberList(channelIndex,sensorIndex), &
326 joTovsPerChannelSensor(channelIndex,sensorIndex)
327 end do loopChannel
328
329 end do loopSensor3
330 write(*,*) ' '
331 end if
332
333 ! print SST data per instrument
334 if(oer_getSSTdataParam_int('numberSSTDatasets') > 0) then
335 write(*,*) 'cfn_sumJo: SST data by data type:'
336 write(*,'(a5,a15,a10,a30,a10, a15)') 'index', ' instrument', ' sensor', ' Jo', ' nobs', ' Jo/nobs'
337 do SSTdatasetIndex = 1, oer_getSSTdataParam_int('numberSSTDatasets')
338 if (nobsInstrumentGlob(SSTdatasetIndex) > 0) then
339 write(*,'(i5,a15,a10,f30.17,i15,f10.5)') SSTdatasetIndex, &
340 oer_getSSTdataParam_char('instrument', SSTdatasetIndex),&
341 oer_getSSTdataParam_char('sensor', SSTdatasetIndex), &
342 joSSTInstrument(SSTdatasetIndex), &
343 nobsInstrumentGlob(SSTdatasetIndex),&
344 joSSTInstrument(SSTdatasetIndex) / &
345 real(nobsInstrumentGlob(SSTdatasetIndex))
346 end if
347 end do
348 end if
349
350 end if
351
352 if(oer_getSSTdataParam_int('numberSSTDatasets') > 0) then
353 deallocate(joSSTInstrument)
354 deallocate(nobsInstrument)
355 deallocate(nobsInstrumentGlob)
356 end if
357
358 end subroutine cfn_sumJo
359
360 !--------------------------------------------------------------------------
361 ! readNameList
362 !--------------------------------------------------------------------------
363 subroutine readNameList
364 !
365 !:Purpose: Reading NAMCFN namelist by any subroutines in costfunction_mod module.
366 !
367 implicit none
368
369 ! Locals:
370 integer :: nulnam, ierr
371 integer, external :: fnom, fclos
372 logical, save :: nmlAlreadyRead = .false.
373 NAMELIST /NAMCFN/ sensorNameList, channelNumberList
374
375 if ( .not. nmlAlreadyRead ) then
376 nmlAlreadyRead = .true.
377
378 !- Setting default values
379 sensorNameList(:) = ''
380 channelNumberList(:,:) = 0
381
382 if ( .not. utl_isNamelistPresent('NAMCFN','./flnml') ) then
383 if ( mmpi_myid == 0 ) then
384 write(*,*) 'NAMCFN is missing in the namelist. The default values will be taken.'
385 end if
386
387 else
388 ! Reading the namelist
389 nulnam = 0
390 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
391 read(nulnam, nml=namcfn, iostat=ierr)
392 if ( ierr /= 0) call utl_abort('costfunction_mod: Error reading namelist')
393 ierr = fclos(nulnam)
394
395 call sortChannelNumbersInNml
396 end if
397 if ( mmpi_myid == 0 ) write(*,nml=namcfn)
398 end if
399
400 end subroutine readNameList
401
402 !--------------------------------------------------------------------------
403 ! sortChannelNumbersInNml
404 !--------------------------------------------------------------------------
405 subroutine sortChannelNumbersInNml
406 !
407 !:Purpose: Sort channelNumbers in NAMCFN namelist. This involves removing
408 ! the duplicates and combine channelNumbers of same sensor
409 ! prescribed on different lines.
410 !
411 implicit none
412
413 ! Locals:
414 integer :: channelIndex
415 integer :: channelIndex1, channelIndex2
416 integer :: sensorIndexInList, sensorIndexInList1, sensorIndexInList2
417
418 character(len=15) :: sensorName1LowerCase, sensorName2LowerCase
419
420 if ( mmpi_myid == 0 ) then
421 write(*,*) 'costfunction_mod: sortChannelNumbersInNml START'
422 end if
423
424 ! set duplicate channelNumber for each sensor to zero
425 loopSensor3: do sensorIndexInList = 1, tvs_nsensors
426 if ( trim(sensorNameList(sensorIndexInList)) == '' ) cycle loopSensor3
427
428 loopChannel2: do channelIndex = tvs_maxNumberOfChannels, 2, -1
429 if ( channelNumberList(channelIndex,sensorIndexInList) == 0 ) cycle loopChannel2
430
431 if ( any(channelNumberList(1:channelIndex-1,sensorIndexInList) == &
432 channelNumberList(channelIndex ,sensorIndexInList)) ) then
433 channelNumberList(channelIndex,sensorIndexInList) = 0
434 end if
435 end do loopChannel2
436 end do loopSensor3
437
438 ! remove later duplicates of sensor/channel pairs
439 loopSensor4: do sensorIndexInList1 = 1, tvs_nsensors-1
440 if ( trim(sensorNameList(sensorIndexInList1)) == '' ) cycle loopSensor4
441
442 call up2low(sensorNameList(sensorIndexInList1),sensorName1LowerCase)
443
444 loopSensor5: do sensorIndexInList2 = sensorIndexInList1+1, tvs_nsensors
445 if ( trim(sensorNameList(sensorIndexInList2)) == '' ) cycle loopSensor5
446
447 call up2low(sensorNameList(sensorIndexInList2),sensorName2LowerCase)
448
449 if ( trim(sensorName2LowerCase) == trim(sensorName1LowerCase) ) then
450 loopChannel3: do channelIndex2 = 1, tvs_maxNumberOfChannels
451 if ( channelNumberList(channelIndex2,sensorIndexInList2) == 0 ) cycle loopChannel3
452
453 if ( any(channelNumberList(: ,sensorIndexInList1) == &
454 channelNumberList(channelIndex2,sensorIndexInList2)) ) then
455 ! set second occurance of channelNumber to zero
456 channelNumberList(channelIndex2,sensorIndexInList2) = 0
457 else
458
459 ! replace the first zero channelNumber of first sensor with new non-zero channelNumber
460 do channelIndex1 = 1, tvs_maxNumberOfChannels
461 if ( channelNumberList(channelIndex1,sensorIndexInList1) == 0 ) then
462 channelNumberList(channelIndex1,sensorIndexInList1) = &
463 channelNumberList(channelIndex2,sensorIndexInList2)
464 channelNumberList(channelIndex2,sensorIndexInList2) = 0
465 exit
466 end if
467 end do
468
469 end if
470
471 end do loopChannel3
472 end if
473 end do loopSensor5
474 end do loopSensor4
475
476 ! if all entries for a sensor are zero, remove the sensor from namelist.
477 ! otherwise sort in ascending order
478 loopSensor6: do sensorIndexInList = 1, tvs_nsensors
479 if ( trim(sensorNameList(sensorIndexInList)) == '' ) cycle loopSensor6
480
481 if ( all(channelNumberList(:,sensorIndexInList) == 0) ) then
482 sensorNameList(sensorIndexInList) = ''
483 else
484 call isort(channelNumberList(:,sensorIndexInList),tvs_maxNumberOfChannels)
485 end if
486 end do loopSensor6
487
488 if ( mmpi_myid == 0 ) then
489 write(*,*) 'costfunction_mod: sortChannelNumbersInNml END'
490 end if
491
492 end subroutine sortChannelNumbersInNml
493
494end module costFunction_mod
