1module obsOperators_mod
2 ! MODULE obsOperators_mod (prefix='oop' category='5. Observation operators')
3 !
4 !:Purpose: All observation operators, including nonlinear, tangent-linear
5 ! and adjoint versions.
6 !
7 use codePrecision_mod
8 use earthConstants_mod
9 use mathPhysConstants_mod
10 use obsSpaceData_mod
11 use columnData_mod
12 use bufr_mod
13 use physicsFunctions_mod
14 use gps_mod
15 use midasMpi_mod
16 use timeCoord_mod
17 use tovsNL_mod
18 use utilities_mod
19 use tovsLin_mod
20 use verticalCoord_mod
21 use varNameList_mod
22 use obsOperatorsChem_mod
23 use obserrors_mod
24 implicit none
25 save
26 private
27
28 ! public procedures
29 public :: oop_ppp_nl, oop_sfc_nl, oop_zzz_nl, oop_gpsro_nl, oop_hydro_nl
30 public :: oop_gpsgb_nl, oop_tovs_nl, oop_chm_nl, oop_sst_nl, oop_ice_nl, oop_raDvel_nl
31 public :: oop_Htl, oop_Had, oop_vobslyrs, oop_iceScaling
32
33 integer, external :: get_max_rss
34
35 real(8), parameter :: temperatureLapseRate = 0.0065D0 ! K/m (i.e. 6.5 K/km)
36
37 ! Jacobian caches
38 real*8 , allocatable :: oop_vRO_Jacobian(:,:,:)
39 logical, allocatable :: oop_vRO_lJac(:)
40 real*8 , allocatable :: oop_vZTD_Jacobian(:,:)
41
42contains
43
44 !--------------------------------------------------------------------------
45 ! oop_vobslyrs
46 !--------------------------------------------------------------------------
47 subroutine oop_vobslyrs( columnTrl, obsSpaceData, beSilent )
48 !:Purpose:
49 ! Find which model levels to use for the vertical interpolation
50 ! of model fields to obs data.
51 implicit none
52
53 ! Arguments:
54 type(struct_columnData), intent(in) :: columnTrl
55 type(struct_obs), intent(inout) :: obsSpaceData
56 logical, intent(in) :: beSilent
57
58 ! Locals:
59 integer :: levIndex,JDATA,NLEV
60 real(8) :: ZLEV,ZPT,ZPB
61 integer :: IOBS,layerIndex,bufrCode
62 character(len=4) :: varLevel
63 integer :: bodyIndex
64
65 if (.not.beSilent) write(*,*) "Entering subroutine OOP_VOBSLYRS"
66
67 ! 2D mode patch
68 if ( col_getNumLev(columnTrl,'MM') <= 1 ) then
69 do bodyIndex = 1, obs_numbody( obsSpaceData )
70 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated .and. &
71 obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) == 1 ) then
72 call obs_bodySet_i(obsSpaceData,OBS_LYR,bodyIndex, 0) ! set OBS_LYR = 0
73 end if
74 end do
75 return
76 end if
77
78 !
79 !-----------------------------------------------------------------------
80 ! --------
81 ! ETA
82 ! --------
83 !
84 ! 1. Find where extrapolation is needed
85 ! ----------------------------------
86 !
87 ! 1.1 PPP Vertical coordinate
88 !
89
90 !$OMP PARALLEL DO PRIVATE(jdata,zlev,iobs,bufrCode,varLevel,zpt,zpb)
91 do JDATA= 1,obs_numbody(obsSpaceData)
92 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,JDATA) == obs_assimilated .and. &
93 obs_bodyElem_i(obsSpaceData,OBS_VCO,JDATA) == 2 ) THEN
94 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,JDATA) /= bufr_nedz ) THEN
95 ZLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,JDATA)
96 else
97 call utl_abort('oop_vobslyr: ZLEV cannot be set, bufr_nedz not supported!')
98 end if
99 IOBS = obs_bodyElem_i(obsSpaceData,OBS_HIND,JDATA)
100 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,JDATA)
101 if (bufr_IsAtmosConstituent(bufrCode)) then
102 varLevel = vnl_varLevelFromVarnum(bufrCode, &
103 obs_headElem_i(obsSpaceData,OBS_CHM,IOBS))
104 else
105 varLevel = vnl_varLevelFromVarnum(bufrCode)
106 end if
107 ZPT= col_getPressure(columnTrl,1,IOBS,varLevel)
108 ZPB= col_getPressure(columnTrl,COL_GETNUMLEV(columnTrl,varLevel),IOBS,varLevel)
109 if ( ZLEV < ZPT ) THEN
110 call obs_bodySet_i(obsSpaceData,OBS_XTR,JDATA,1)
111 !
112 !- !!! WARNING !!! This obs is above the model lid.
113 ! We must turn off its assimilation flag because the
114 ! current obs operators cannot deal with this situation (JFC)
115 if (varLevel /= 'SF') then
116 write(*,*) 'oop_vobslyrs: Rejecting OBS above model lid, pressure = ', ZLEV,' < ',ZPT
117 call obs_bodySet_i(obsSpaceData,OBS_ASS,JDATA, obs_notAssimilated)
118 end if
119 else if ( ZLEV > ZPB ) THEN
120 call obs_bodySet_i(obsSpaceData,OBS_XTR,JDATA,2)
121 else
122 call obs_bodySet_i(obsSpaceData,OBS_XTR,JDATA,0)
123 end if
124 end if
125 end do
126 !$OMP END PARALLEL DO
127 !
128 ! 1.2 ZZZ Vertical coordinate
129 !
130 !$OMP PARALLEL do PRIVATE(jdata,zlev,iobs,bufrCode,varLevel,zpt,zpb,nlev)
131 do JDATA= 1,obs_numbody(obsSpaceData)
132 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,JDATA) == obs_assimilated .and. &
133 obs_bodyElem_i(obsSpaceData,OBS_VCO,JDATA) == 1 ) then
134 IOBS = obs_bodyElem_i(obsSpaceData,OBS_HIND,JDATA)
135 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,JDATA)
136 if ( bufrCode /= bufr_nedz ) then
137 ZLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,JDATA)
138 if ( bufrCode == bufr_nebd ) THEN
139 ZLEV = ZLEV - obs_headElem_r(obsSpaceData,OBS_TRAD,IOBS)
140 endif
141 else
142 call utl_abort('oop_vobslyr: ZLEV cannot be set, bufr_nedz not supported!')
143 end if
144 if (bufr_IsAtmosConstituent(bufrCode)) then
145 varLevel = vnl_varLevelFromVarnum(bufrCode, &
146 obs_headElem_i(obsSpaceData,OBS_CHM,IOBS))
147 else
148 varLevel = vnl_varLevelFromVarnum(bufrCode)
149 end if
150 if (varLevel == 'SF') then
151 ZPT= col_getHeight(columnTrl,1,IOBS,'TH')
152 ZPB= col_getHeight(columnTrl,0,IOBS,'SF')
153 else
154 nlev=col_getNumLev(columnTrl,varLevel)
155 ZPT= col_getHeight(columnTrl,1,IOBS,varLevel)
156 ZPB= col_getHeight(columnTrl,NLEV,IOBS,varLevel)
157 end if
158 if ( ZLEV > ZPT ) then
159 call obs_bodySet_i(obsSpaceData,OBS_XTR,JDATA,1)
160 write(*,*) 'oop_vobslyrs: Rejecting OBS above model lid, height =', ZLEV,' > ',ZPT
161 call obs_bodySet_i(obsSpaceData,OBS_ASS,JDATA, obs_notAssimilated)
162 else if ( ZLEV < ZPB ) then
163 call obs_bodySet_i(obsSpaceData,OBS_XTR,JDATA,2)
164 else
165 call obs_bodySet_i(obsSpaceData,OBS_XTR,JDATA,0)
166 end if
167 end if
168 end do
169 !$OMP END PARALLEL DO
170 !
171 !
172 ! 2. FInd interpolation layer
173 ! ------------------------
174 ! (Model levels are assumed to be in increasing order in Mbs)
175 !
176 ! 2.1 PPP Vertical coordinate
177 !
178 !$OMP PARALLEL DO PRIVATE(jdata,iobs,zlev,bufrCode,varLevel,layerIndex,nlev,levIndex,zpt,zpb)
179 do JDATA = 1, obs_numbody(obsSpaceData)
180 call obs_bodySet_i(obsSpaceData,OBS_LYR,JDATA,0)
181 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,JDATA) == obs_assimilated .and. &
182 obs_bodyElem_i(obsSpaceData,OBS_VCO,JDATA) == 2 ) then
183 IOBS = obs_bodyElem_i(obsSpaceData,OBS_HIND,JDATA)
184 ZLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,JDATA)
185 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,JDATA)
186 if (bufr_IsAtmosConstituent(bufrCode)) then
187 varLevel = vnl_varLevelFromVarnum(bufrCode, &
188 obs_headElem_i(obsSpaceData,OBS_CHM,IOBS))
189 else
190 varLevel = vnl_varLevelFromVarnum(bufrCode)
191 end if
192 layerIndex = 1
193 nlev=COL_GETNUMLEV(columnTrl,varLevel)
194 do levIndex = 2,NLEV - 1
195 ZPT = col_getPressure(columnTrl,levIndex,IOBS,varLevel)
196 if ( ZLEV > ZPT ) layerIndex = levIndex
197 end do
198 ZPT = col_getPressure(columnTrl,layerIndex,IOBS,varLevel)
199 ZPB = col_getPressure(columnTrl,layerIndex+1,IOBS,varLevel)
200 call obs_bodySet_i(obsSpaceData,OBS_LYR,JDATA, layerIndex)
201 end if
202 end do
203 !$OMP END PARALLEL DO
204 !
205 ! 2.2 ZZZ Vertical coordinate and surface observations
206 !
207 !$OMP PARALLEL DO PRIVATE(jdata,iobs,zlev,bufrCode,varLevel,layerIndex,nlev,levIndex,zpt)
208 do JDATA = 1, obs_numbody(obsSpaceData)
209 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,JDATA) == obs_assimilated .and. &
210 obs_bodyElem_i(obsSpaceData,OBS_VCO,JDATA) == 1 ) then
211 IOBS = obs_bodyElem_i(obsSpaceData,OBS_HIND,JDATA)
212 ZLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,JDATA)
213 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,JDATA)
214 if (bufr_IsAtmosConstituent(bufrCode)) then
215 varLevel = vnl_varLevelFromVarnum(bufrCode, &
216 obs_headElem_i(obsSpaceData,OBS_CHM,IOBS))
217 else
218 varLevel = vnl_varLevelFromVarnum(bufrCode)
219 end if
220 layerIndex = 1
221 nlev=COL_GETNUMLEV(columnTrl,varLevel)
222 do levIndex = 2, NLEV - 1
223 ZPT = col_getHeight(columnTrl,levIndex,IOBS,varLevel)
224 if ( ZLEV < ZPT ) layerIndex = levIndex
225 end do
226 if ( bufrCode == bufr_neps .or. bufrCode == bufr_nepn .or. &
227 bufrCode == bufr_nezd .or. bufrCode == bufr_gust .or. &
228 bufrCode == bufr_radarPrecip .or. bufrCode == bufr_logRadarPrecip ) THEN
229 ! for surface observations associated with surface analysis variables
230 layerIndex = 0
231 else if ( bufrCode == bufr_nets .or. bufrCode == bufr_ness .or. &
232 bufrCode == bufr_neus .or. bufrCode == bufr_nevs .or. &
233 bufrCode == bufr_nehs .or. bufrCode == bufr_vis .or. &
234 bufrCode == bufr_logVis ) then
235 ! for surface observations associated with NON-surface analysis variables
236 layerIndex = nlev - 1
237 end if
238 call obs_bodySet_i(obsSpaceData,OBS_LYR,JDATA, layerIndex)
239 end if
240 end do
241 !$OMP END PARALLEL DO
242 !
243 end subroutine oop_vobslyrs
244
245 !--------------------------------------------------------------------------
246 ! oop_ppp_nl
247 !--------------------------------------------------------------------------
248 subroutine oop_ppp_nl( columnTrlOnTrlLev, obsSpaceData, beSilent, cdfam, destObsColumn )
249 !:Purpose: Computation of y - H(x)
250 ! for pressure-level observations.
251 ! Interpolate vertically columnTrlOnTrlLev to
252 ! the pressure levels of the observations.
253 ! A linear interpolation in ln(p) is performed.
254 !
255 implicit none
256
257 ! Arguments:
258 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
259 type(struct_obs) , intent(inout) :: obsSpaceData
260 logical , intent(in) :: beSilent
261 character(len=*) , intent(in) :: cdfam ! family of obsservation
262 integer , intent(in) :: destObsColumn
263
264 ! Locals:
265 integer :: headerIndex,bodyIndex,ilyr
266 integer :: iass,ixtr,ivco,bufrCode,nlev_T
267 real(8) :: zvar,zwb,zwt,zexp
268 real(8) :: zlev,zpt,zpb,zomp,ztvg
269 real(8) :: trlValueBot,trlValueTop,lat
270 character(len=4) :: varName
271 character(len=4) :: varLevel
272 real(8),pointer :: col_ptr(:),col_ptr_tt(:),col_ptr_hu(:)
273 real(8), allocatable :: geopotential(:)
274 real(8) :: heightSfc(1), geopotentialSfc(1)
275
276 if (.not.beSilent) write(*,*) 'Entering subroutine oop_ppp_nl'
277
278 zexp = MPC_RGAS_DRY_AIR_R8 * temperatureLapseRate / ec_rg
279
280 nlev_T = col_getNumLev(columnTrlOnTrlLev,'TH')
281 allocate(geopotential(nlev_T))
282
283 call obs_set_current_body_list(obsSpaceData, cdfam)
284 BODY: do
285 bodyIndex = obs_getBodyIndex(obsSpaceData)
286 if (bodyIndex < 0) exit BODY
287
288 ! Only process pressure level observations flagged to be assimilated
289 iass=obs_bodyElem_i (obsSpaceData,OBS_ASS,bodyIndex)
290 ivco=obs_bodyElem_i (obsSpaceData,OBS_VCO,bodyIndex)
291 if (iass /= 1 .or. ivco /= 2) cycle BODY
292
293 ixtr=obs_bodyElem_i (obsSpaceData,OBS_XTR,bodyIndex)
294 bufrCode=obs_bodyElem_i (obsSpaceData,OBS_VNM,bodyIndex)
295 zvar=obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
296 zlev=obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
297 headerIndex=obs_bodyElem_i (obsSpaceData,OBS_HIND,bodyIndex)
298
299 if ( ixtr == 0 ) then
300
301 ! Process all data within the domain of the model
302 ilyr =obs_bodyElem_i (obsSpaceData,OBS_LYR,bodyIndex)
303 varName = vnl_varNameFromVarnum(bufrCode)
304 varLevel = vnl_varLevelFromVarnum(bufrCode)
305 zpt= col_getPressure(columnTrlOnTrlLev,ilyr ,headerIndex,varLevel)
306 zpb= col_getPressure(columnTrlOnTrlLev,ilyr+1,headerIndex,varLevel)
307 zwb = log(zlev/zpt)/log(zpb/zpt)
308 zwt = 1.d0 - zwb
309 lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
310 if (bufrCode == bufr_nees) then
311 col_ptr_hu=>col_getColumn(columnTrlOnTrlLev,headerIndex,'HU')
312 col_ptr_tt=>col_getColumn(columnTrlOnTrlLev,headerIndex,'TT')
313 trlValueBot=hutoes(col_ptr_hu(ilyr+1),col_ptr_tt(ilyr+1),zpb)
314 trlValueTop=hutoes(col_ptr_hu(ilyr ),col_ptr_tt(ilyr ),zpt)
315 else
316 if (trim(varName) == 'Z_T') then
317 col_ptr=>col_getColumn(columnTrlOnTrlLev,headerIndex,'Z_T')
318 call phf_height2geopotential(col_ptr,lat,geopotential)
319 trlValueBot=geopotential(ilyr+1)
320 trlValueTop=geopotential(ilyr )
321 else
322 col_ptr=>col_getColumn(columnTrlOnTrlLev,headerIndex,varName)
323 trlValueBot=col_ptr(ilyr+1)
324 trlValueTop=col_ptr(ilyr )
325 end if
326 end if
327 zomp = zvar-(zwb*trlValueBot+zwt*trlValueTop)
328 call obs_bodySet_r(obsSpaceData,destObsColumn,bodyIndex,zomp)
329
330 else if (ixtr == 2) then
331
332 ! Process only GZ that is data below model's orography
333 if (bufrCode == bufr_negz ) then
334 ! Forward nonlinear model for geopotential data below model's orography
335 ztvg = (1.0d0 + MPC_DELTA_R8 * col_getElem(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'TH'),headerIndex,'HU'))* &
336 col_getElem(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'TH'),headerIndex,'TT')
337
338 ! convert height of surface to geopotential
339 heightSfc(1) = col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF')
340 lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
341 call phf_height2geopotential(heightSfc,lat,geopotentialSfc)
342
343 zomp = ( zvar - geopotentialSfc(1) - &
344 ztvg/(temperatureLapseRate/ec_rg)*(1.D0-(zlev/col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0'))**zexp))
345 call obs_bodySet_r(obsSpaceData,destObsColumn,bodyIndex,zomp)
346 end if
347
348 end if
349
350 end do body
351
352 deallocate(geopotential)
353
354 end subroutine oop_ppp_nl
355
356 !--------------------------------------------------------------------------
357 ! oop_zzz_nl
358 !--------------------------------------------------------------------------
359 subroutine oop_zzz_nl( columnTrlOnTrlLev, obsSpaceData, beSilent, cdfam, &
360 destObsColumn )
361 !:Purpose: Computation of y - H(x) for geometric-height observations
362 ! Interpolate vertically columnTrlOnTrlLev to the geometric heights (in
363 ! meters) of the observations.
364 ! A linear interpolation in z is performed.
365 !
366 !:Notes:
367 ! As a first approximation, use the geopotential height. Once this is
368 ! working, this should be changed for a calculation of the geometric
369 ! height.
370 !
371 ! Note that, in the case of an aladin HLOS wind, the correction to zvar
372 ! (OBS_VAR) is not written back to obsSpaceData. It is simply used to
373 ! calculate OMP (which is written to obsSpaceData) and then is discarded.
374 ! Thereafter, if one calculates OMP - O (this will be the uncorrected O),
375 ! the result will be a corrected P.
376 implicit none
377
378 ! Arguments:
379 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
380 type(struct_obs), intent(inout) :: obsSpaceData
381 logical, intent(in) :: beSilent
382 integer, intent(in) :: destObsColumn
383 character(len=*), intent(in) :: cdfam ! family of observation
384
385 ! Locals:
386 integer :: headerIndex,bodyIndex,ilyr,bufrCode,levIndexTop,levIndexBot
387 integer :: bodyIndexStart,bodyIndexEnd,bodyIndex2
388 integer :: found ! a group of bit flags
389 integer :: ierr, nulnam, fnom,fclos
390 real(8) :: zvar,zwb,zwt
391 real(8) :: zlev,zpt,zpb,zomp
392 real(8) :: trlValueBot,trlValueTop
393 character(len=4) :: varLevel
394 real(8) :: value ! temporary holder
395 real(8) :: azimuth ! HLOS wind direction CW from true north
396 real(8) :: tempRef ! reference temperature used to calculate HLOS wind
397 real(8) :: presRef ! reference pressure used to calculate HLOS wind
398 real(8) :: dwdp ! derivative of HLOS wind wrt P
399 real(8) :: dwdt ! derivative of HLOS wind wrt T
400 real(8) :: uuLyr, vvLyr ! wind on layer, OBS_LYR
401 real(8) :: uuLyr1,vvLyr1 ! wind on layer plus 1
402 real(8) :: ttLyr, ppLyr ! T, P on layer, OBS_LYR
403 real(8) :: ttLyr1,ppLyr1 ! T, P on layer, OBS_LYR plus 1
404 real(8) :: ttbg, ppbg ! background T, P at the observation location
405 logical :: list_is_empty
406
407 ! Namelist variables:
408 logical :: do_adjust_aladin ! choose to adjust obs value as if it was retrieved using our temp and pressure
409
410 namelist /NAMALADIN_OBS/do_adjust_aladin
411
412 if (.not.beSilent) write(*,*) "Entering subroutine oop_zzz_nl"
413
414 call obs_set_current_body_list(obsSpaceData, cdfam, list_is_empty)
415
416 if (list_is_empty)then
417 return
418 end if
419
420 ! Read in the namelist NAMALADIN_OBS
421 do_adjust_aladin = .false.
422 nulnam=0
423 ierr=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
424 read(nulnam,nml=namaladin_obs,iostat=ierr)
425 if (ierr.ne.0) call utl_abort('oop_zzz_nl: Error reading namelist')
426 if (.not.beSilent) write(*,nml=namaladin_obs)
427 ierr=fclos(nulnam)
428
429 BODY: do
430 bodyIndex = obs_getBodyIndex(obsSpaceData)
431 if (bodyIndex < 0) exit BODY
432
433 ! Process all geometric-height data within the domain of the model
434 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) /= obs_assimilated .or. &
435 obs_bodyElem_i(obsSpaceData,OBS_XTR,bodyIndex) /= 0 .or. &
436 obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) /= 1 ) &
437 cycle BODY
438 ! So, OBS_VCO==1 => OBS_PPP is a height in m
439
440 bufrCode=obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
441 zvar=obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
442 zlev=obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
443 headerIndex=obs_bodyElem_i(obsSpaceData,OBS_HIND,bodyIndex)
444
445 ilyr = obs_bodyElem_i(obsSpaceData,OBS_LYR,bodyIndex)
446 varLevel = vnl_varLevelFromVarnum(bufrCode)
447 zpt= col_getHeight(columnTrlOnTrlLev,ilyr ,headerIndex,varLevel)
448 zpb= col_getHeight(columnTrlOnTrlLev,ilyr+1,headerIndex,varLevel)
449 zwb = (zpt-zlev)/(zpt-zpb)
450 zwt = 1.d0 - zwb
451
452 select case (bufrCode)
453 case (bufr_neal) ! Aladin HLOS wind observation
454 ! Scan body indices for the needed attributes
455 found = 0
456 bodyIndexStart = obs_headElem_i(obsSpaceData, OBS_RLN, headerIndex)
457 bodyIndexEnd = obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex) &
458 + bodyIndexStart - 1
459 tempRef = 0.0d0
460 dwdt = 0.0d0
461 presRef = 0.0d0
462 dwdp = 0.0d0
463 BODY_SUPP: do bodyIndex2 = bodyIndexStart, bodyIndexEnd
464
465 value = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex2)
466 select case(obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex2))
467 case(bufr_NEAZ)
468 azimuth = value * MPC_RADIANS_PER_DEGREE_R8
469 found = ibset(found,0)
470
471 case(bufr_nett)
472 tempRef = value
473 found = ibset(found,1)
474
475 case(bufr_neps)
476 presRef = value
477 found = ibset(found,2)
478
479 case(bufr_NEDWDP)
480 dwdp = value
481 found = ibset(found,3)
482
483 case(bufr_NEDWDT)
484 dwdt = value
485 found = ibset(found,4)
486 end select
487
488 if (popcnt(found) == 5) exit BODY_SUPP
489 end do BODY_SUPP
490
491 if (.not. btest(found,0))then
492 ! The azimuth was not found. The observation cannot be treated
493 ! Set the assimilation flag to 0 to ignore this datum later.
494 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
495 cycle BODY
496 end if
497
498 ! Obtain the needed forecast data
499 uuLyr =col_getElem(columnTrlOnTrlLev,ilyr, headerIndex,'UU')
500 uuLyr1=col_getElem(columnTrlOnTrlLev,ilyr+1,headerIndex,'UU')
501 vvLyr =col_getElem(columnTrlOnTrlLev,ilyr, headerIndex,'VV')
502 vvLyr1=col_getElem(columnTrlOnTrlLev,ilyr+1,headerIndex,'VV')
503 ttLyr =col_getElem(columnTrlOnTrlLev,ilyr, headerIndex,'TT')
504 ttLyr1=col_getElem(columnTrlOnTrlLev,ilyr+1,headerIndex,'TT')
505 ppLyr =col_getPressure(columnTrlOnTrlLev,ilyr ,headerIndex,'MM')
506 ppLyr1=col_getPressure(columnTrlOnTrlLev,ilyr+1,headerIndex,'MM')
507
508 ! Interpolate forecast T, P to the observation location
509 ttbg = zwb*ttLyr1 + zwt*ttLyr
510 ppbg = zwb*ppLyr1 + zwt*ppLyr
511
512 ! Adjust zvar, the HLOS wind observation, if all attributes are available
513 if (do_adjust_aladin .and. (popcnt(found) == 5)) then
514 ! Adjust in situ the HLOS wind data from obsSpaceData to account for
515 ! the differences between our T, P forecast fields and those of the NWP
516 ! site that calculated the HLOS wind values. The goal is to produce an
517 ! HLOS wind observation as if it had been calculated by us.
518 zvar = zvar + (ttbg - tempRef) * dwdt &
519 + (ppbg - presRef) * dwdp
520 end if
521
522 ! Apply the nonlinear aladin observation operator
523 trlValueBot= -vvLyr1*cos(azimuth) - uuLyr1*sin(azimuth)
524 trlValueTop= -vvLyr *cos(azimuth) - uuLyr *sin(azimuth)
525
526 ! For aladin data, the temperature and pressure are really *reference*
527 ! values. They must not be assimilated. Mark them so.
528 case (bufr_nett)
529 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
530 cycle BODY
531 case (bufr_neps)
532 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
533 cycle BODY
534 case (bufr_nees)
535 call utl_abort('oop_zzz_nl: CANNOT ASSIMILATE ES!!!')
536
537 case default
538 ! These are the profiler observations
539 levIndexTop = ilyr + col_getOffsetFromVarno(columnTrlOnTrlLev,bufrCode)
540 levIndexBot = levIndexTop+1
541 trlValueBot=col_getElem(columnTrlOnTrlLev,levIndexBot,headerIndex)
542 trlValueTop=col_getElem(columnTrlOnTrlLev,levIndexTop,headerIndex)
543 end select
544
545 zomp = zvar-(zwb*trlValueBot+zwt*trlValueTop)
546 call obs_bodySet_r(obsSpaceData,destObsColumn,bodyIndex,zomp)
547
548 enddo BODY
549
550 end subroutine oop_zzz_nl
551
552 !--------------------------------------------------------------------------
553 ! oop_sfc_nl
554 !--------------------------------------------------------------------------
555 subroutine oop_sfc_nl( columnTrlOnTrlLev, obsSpaceData, beSilent, cdfam, &
556 destObsColumn )
557 !:Purpose: Computation of the residuals to the observations
558 ! FOR SURFACE DATA (except ground-based GPS zenith delay).
559 implicit none
560
561 ! Arguments:
562 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
563 type(struct_obs) , intent(inout) :: obsSpaceData
564 logical , intent(in) :: beSilent
565 character(len=*) , intent(in) :: cdfam ! family of observation
566 integer , intent(in) :: destObsColumn
567
568 ! Locals:
569 integer :: columnLevelIndex, bufrCode, headerIndex, bodyIndex
570 integer :: ierr, nulnam, fnom, fclos
571 real(8) :: obsValue, trlVirtTemp
572 real(8) :: pCorrectionFactor, coeffA, coeffB
573 real(8) :: obsHeight, deltaT, delTdelZ, trlLevelHeight
574 real(8) :: trlValue
575 real(8) :: trlUwind, trlVwind, squareSum, trlWindSpeed
576 character(len=4) :: varLevel
577 logical, save :: firstCall = .true.
578
579 ! Namelist variables:
580 logical, save :: adjustTemperature ! choose to adjust near-sfc temperature using lapse rate and height difference
581
582 namelist /namSurfaceObs/adjustTemperature
583
584 if (.not.beSilent) write(*,*) "Entering subroutine oop_sfc_nl"
585
586 ! Read in the namelist namSurfaceObs
587 if (firstCall) then
588 adjustTemperature = .true. ! default value
589
590 if (utl_isNamelistPresent('namSurfaceObs','./flnml')) then
591 nulnam=0
592 ierr=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
593 read(nulnam,nml=namSurfaceObs,iostat=ierr)
594 if (ierr /= 0) call utl_abort('oop_sfc_nl: Error reading namelist namSurfaceObs')
595 if (.not. beSilent) write(*,nml=namSurfaceObs)
596 ierr=fclos(nulnam)
597 else if (.not. beSilent) then
598 write(*,*)
599 write(*,*) 'oop_sfc_nl: namSurfaceObs is missing in the namelist. The default value will be taken.'
600 end if
601 firstCall = .false.
602 end if
603
604 ! loop over all header indices of the specified family with surface obs
605 call obs_set_current_header_list(obsSpaceData,cdfam)
606 HEADER: do
607 headerIndex = obs_getHeaderIndex(obsSpaceData)
608 if (headerIndex < 0) exit HEADER
609 ! loop over all body indices for this headerIndex
610 call obs_set_current_body_list(obsSpaceData, headerIndex)
611 BODY: do
612 bodyIndex = obs_getBodyIndex(obsSpaceData)
613 if (bodyIndex < 0) exit BODY
614
615 ! only process height level observations flagged to be assimilated
616 if ( obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) /= 1 .or. &
617 obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) /= obs_assimilated ) cycle BODY
618
619 ! only process this set of surface observations
620 bufrCode=obs_bodyElem_i (obsSpaceData,OBS_VNM,bodyIndex)
621 if (bufrCode /= bufr_nets .and. bufrCode /= bufr_neps .and. &
622 bufrCode /= bufr_neus .and. bufrCode /= bufr_nevs .and. &
623 bufrCode /= bufr_ness .and. bufrCode /= bufr_nepn .and. &
624 bufrCode /= bufr_vis .and. bufrCode /= bufr_logVis .and. &
625 bufrCode /= bufr_gust .and. &
626 bufrCode /= bufr_radarPrecip .and. &
627 bufrCode /= bufr_logRadarPrecip .and. &
628 bufrCode /= bufr_nefs ) cycle BODY
629
630 obsValue = obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
631 obsHeight= obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
632 ! obsHeigth = station height + obs validity height offset as defined in obsUtil_mod->surfvcord
633 varLevel = vnl_varLevelFromVarnum(bufrCode)
634
635 if (bufrCode == bufr_nets .or. bufrCode == bufr_ness .or. &
636 bufrCode == bufr_neus .or. bufrCode == bufr_nevs .or. &
637 bufrCode == bufr_gust .or. bufrCode == bufr_vis .or. &
638 bufrCode == bufr_logVis .or. &
639 bufrCode == bufr_radarPrecip .or. &
640 bufrCode == bufr_logRadarPrecip ) then
641
642 ! T1.5m,(T-Td)1.5m,u10m,v10m
643 ! In this section we always extrapolate linearly the trial
644 ! field at the model surface to the height of the
645 ! surface observation whether the observation is above or
646 ! below the model surface.
647 ! NOTE: For (T-Td)1.5m,u10m,v10m we do a zero order extrapolation
648
649 if (bufrCode == bufr_nets .and. adjustTemperature) then
650 delTdelZ = temperatureLapseRate
651 else
652 delTdelZ = 0.0d0
653 end if
654
655 if (bufrCode == bufr_ness) then
656 trlValue = hutoes(col_getElem(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'TH'),headerIndex,'HU'), &
657 col_getElem(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'TH'),headerIndex,'TT'), &
658 col_getPressure(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'TH'),headerIndex,'TH'))
659 else
660 columnLevelIndex = col_getNumLev(columnTrlOnTrlLev,varLevel) + col_getOffsetFromVarno(columnTrlOnTrlLev,bufrCode)
661 trlValue = col_getElem(columnTrlOnTrlLev,columnLevelIndex,headerIndex)
662 end if
663 trlLevelHeight = col_getHeight(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,varLevel),headerIndex,varLevel)
664
665 call obs_bodySet_r(obsSpaceData,destObsColumn,bodyIndex, &
666 (obsValue-trlValue + delTdelZ*(obsHeight-trlLevelHeight)) )
667
668 else if ( bufrCode == bufr_neps .or. bufrCode == bufr_nepn ) then
669
670 ! Surface (PS) & mean sea level (PN) pressure cases
671 ! Background surface pressure are corrected for the height difference with the
672 ! observation. For mean sea level observation, the observation height = 0.
673
674 ! 1) Temperature difference = lapse-rate (6.5 degree/km) * height difference (dz)
675 deltaT = temperatureLapseRate*(obsHeight-col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF'))
676
677 ! 2) Compute the 1.5m background virtual temperature: Tv = T*(1+0.608*HU)
678 trlVirtTemp = (1.0d0 + MPC_DELTA_R8 * &
679 col_getElem(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'TH'),headerIndex,'HU')) * &
680 col_getElem(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'TH'),headerIndex,'TT')
681
682 ! 3) Compute the correction coefficient
683 ! The legacy code says...
684 coeffA = (trlVirtTemp-deltaT)/trlVirtTemp
685 coeffB = 1.0D0/(MPC_RGAS_DRY_AIR_R8*temperatureLapseRate/ec_rg)
686 pCorrectionFactor = coeffA**coeffB
687 ! However, the U.S. Standard Atmosphere (1976, U.S. Government Printing Office, Washington, D.C., 1976*)
688 ! at page 12 says...
689 ! coeffA = trlVirtTemp/(trlVirtTemp+deltaT)
690 ! but the former was found to perform better (gives lower O-P values) than the latter by J-F Caron in 2018
691
692 ! 4) O-P, where P = P0 * pCorrectionFactor (Eq. 33a at page 12 of the U.S. Standard Atmosphere, 1976,
693 ! U.S. Government Printing Office, Washington, D.C., 1976*)
694 call obs_bodySet_r(obsSpaceData,destObsColumn,bodyIndex, &
695 obsValue-(col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0')*pCorrectionFactor))
696
697 ! (*) available at https://ntrs.nasa.gov/api/citations/19770009539/downloads/19770009539.pdf
698
699 else if (bufrCode == bufr_nefs) then
700
701 trlUwind = col_getElem(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'MM'),headerIndex,'UU')
702 trlVwind = col_getElem(columnTrlOnTrlLev,col_getNumLev(columnTrlOnTrlLev,'MM'),headerIndex,'VV')
703 squareSum = trlUwind**2 + trlVwind**2
704 if ( squareSum > 1.d-10 ) then
705 trlWindSpeed = sqrt(squareSum)
706 else
707 trlWindSpeed = 0.0
708 end if
709 call obs_bodySet_r(obsSpaceData,OBS_OMP,bodyIndex, &
710 obsValue-trlWindSpeed)
711
712 end if
713
714 end do BODY
715
716 end do HEADER
717
718 end subroutine oop_sfc_nl
719
720 !--------------------------------------------------------------------------
721 ! oop_sst_nl
722 !--------------------------------------------------------------------------
723 subroutine oop_sst_nl( columnTrlOnTrlLev, obsSpaceData, beSilent, cdfam, &
724 destObsColumn )
725 !:Purpose: Computation of Jo and the residuals to the observations
726 ! for Sea Surface Temperature (SST) data.
727 implicit none
728
729 ! Arguments:
730 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
731 type(struct_obs) , intent(inout) :: obsSpaceData
732 logical , intent(in) :: beSilent
733 character(len=*) , intent(in) :: cdfam ! family of observation
734 integer , intent(in) :: destObsColumn
735
736 ! Locals:
737 integer :: bufrCode, headerIndex, bodyIndex
738 real(8) :: obsValue
739 character(len=4) :: varName
740
741 if (.not.beSilent) write(*,*) "Entering subroutine oop_sst_nl, family: ", trim(cdfam)
742
743 ! loop over all header indices of the specified family with surface obs
744 call obs_set_current_header_list( obsSpaceData, cdfam )
745
746 HEADER: do
747
748 headerIndex = obs_getHeaderIndex( obsSpaceData )
749 if ( headerIndex < 0 ) exit HEADER
750
751 ! loop over all body indices for this headerIndex
752 call obs_set_current_body_list( obsSpaceData, headerIndex )
753
754 BODY: do
755
756 bodyIndex = obs_getBodyIndex( obsSpaceData )
757 if ( bodyIndex < 0 ) exit BODY
758
759 ! only process observations flagged to be assimilated
760 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) /= obs_assimilated ) cycle BODY
761
762 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
763
764 if ( bufrCode /= bufr_sst ) cycle BODY
765
766 if ( col_varExist(columnTrlOnTrlLev,'TM') ) then
767 varName = 'TM'
768 else
769 varName = 'TG'
770 end if
771
772 obsValue = obs_bodyElem_r( obsSpaceData, OBS_VAR, bodyIndex )
773 call obs_bodySet_r( obsSpaceData, destObsColumn, bodyIndex, &
774 obsValue - ( col_getElem( columnTrlOnTrlLev, 1, headerIndex, varName ) ))
775
776 end do BODY
777
778 end do HEADER
779
780 end subroutine oop_sst_nl
781
782 !--------------------------------------------------------------------------
783 ! oop_hydro_nl
784 !--------------------------------------------------------------------------
785 subroutine oop_hydro_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, cdfam, &
786 destObsColumn)
787 !:Purpose: To computate Jo and the residuals to the observations
788 ! for hydrological data
789 implicit none
790
791 ! Arguments:
792 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
793 type(struct_obs) , intent(inout) :: obsSpaceData
794 logical , intent(in) :: beSilent
795 character(len=*) , intent(in) :: cdfam ! family of observation
796 integer , intent(in) :: destObsColumn
797
798 ! Locals:
799 integer :: bufrCode, headerIndex, bodyIndex
800 real(8) :: obsValue
801 character(len=4) :: varName
802
803 if (.not.beSilent) write(*,*) "Entering subroutine oop_hydro_nl, family: ", trim(cdfam)
804
805 ! loop over all header indices of the specified family with surface obs
806 call obs_set_current_header_list( obsSpaceData, cdfam )
807
808 HEADER: do
809 headerIndex = obs_getHeaderIndex( obsSpaceData )
810 if ( headerIndex < 0 ) exit HEADER
811
812 ! loop over all body indices for this headerIndex
813 call obs_set_current_body_list( obsSpaceData, headerIndex )
814
815 BODY: do
816 bodyIndex = obs_getBodyIndex( obsSpaceData )
817 if ( bodyIndex < 0 ) exit BODY
818
819 ! only process observations flagged to be assimilated
820 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) /= obs_assimilated ) cycle BODY
821
822 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
823
824 if ( bufrCode /= bufr_riverFlow ) cycle BODY
825
826 obsValue = obs_bodyElem_r( obsSpaceData, OBS_VAR, bodyIndex )
827 varName = vnl_varNameFromVarNum(bufrCode)
828 call obs_bodySet_r( obsSpaceData, destObsColumn, bodyIndex, &
829 obsValue - col_getElem(columnTrlOnTrlLev,1,headerIndex, varName_opt = varName) )
830
831 end do BODY
832
833 end do HEADER
834
835 end subroutine oop_hydro_nl
836
837 !--------------------------------------------------------------------------
838 ! oop_ice_nl
839 !--------------------------------------------------------------------------
840 subroutine oop_ice_nl( columnTrlOnTrlLev, obsSpaceData, beSilent, cdfam, &
841 destObsColumn )
842 !:Purpose: Computation of Jo and the residuals to the observations
843 ! FOR SEA ICE CONCENTRATION DATA
844 implicit none
845
846 ! Arguments:
847 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
848 type(struct_obs) , intent(inout) :: obsSpaceData
849 logical , intent(in) :: beSilent
850 character(len=*) , intent(in) :: cdfam ! family of observation
851 integer , intent(in) :: destObsColumn
852
853 ! Locals:
854 integer :: bufrCode, headerIndex, bodyIndex
855 integer :: obsDate, monthIndex
856 integer :: trackCellNum
857 real(8) :: obsValue, backValue
858 real(8) :: conc
859 character(len=4) :: varName
860 character(len=8) :: ccyymmdd
861
862 if (.not. beSilent) write(*,*) "Entering subroutine oop_ice_nl, family: ", trim(cdfam)
863
864 ! loop over all body indices
865 call obs_set_current_body_list( obsSpaceData, cdfam )
866
867 BODY: do
868
869 bodyIndex = obs_getBodyIndex( obsSpaceData )
870 if ( bodyIndex < 0 ) exit BODY
871
872 ! only process observations flagged to be assimilated
873 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) /= obs_assimilated ) cycle BODY
874
875 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
876
877 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
878 varName = vnl_varNameFromVarNum(bufrCode)
879
880 select case (bufrCode)
881 case(bufr_icec, bufr_icep)
882 backValue = 100.0d0*col_getElem( columnTrlOnTrlLev, 1, headerIndex, varName )
883 case(bufr_icev)
884 backValue = 1.0d0*col_getElem( columnTrlOnTrlLev, 1, headerIndex, varName )
885 case(bufr_ices)
886 obsDate = obs_headElem_i( obsSpaceData, OBS_DAT, headerIndex )
887 write(ccyymmdd, FMT='(i8.8)') obsDate
888 read(ccyymmdd(5:6), FMT='(i2)') monthIndex
889 conc = col_getElem( columnTrlOnTrlLev, 1, headerIndex, varName)
890 trackCellNum = obs_headElem_i( obsSpaceData, OBS_FOV, headerIndex )
891 backValue = (1.0d0-conc)*oer_ascatAnisOpenWater(trackCellNum,monthIndex) + &
892 conc*oer_ascatAnisIce(trackCellNum,monthIndex)
893 case default
894 cycle BODY
895 end select
896
897 obsValue = obs_bodyElem_r( obsSpaceData, OBS_VAR, bodyIndex )
898 call obs_bodySet_r( obsSpaceData, destObsColumn, bodyIndex, obsValue - backValue )
899
900 end do BODY
901
902 end subroutine oop_ice_nl
903
904 !--------------------------------------------------------------------------
905 ! oop_raDvel_nl
906 !--------------------------------------------------------------------------
907 subroutine oop_raDvel_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, cdfam, &
908 destObsColumn)
909 !:Purpose: Computation of Jo and OMP for Radar Doppler velocity observations
910 implicit none
911
912 ! Arguments:
913 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
914 type(struct_obs) , intent(inout) :: obsSpaceData
915 logical , intent(in) :: beSilent
916 character(len=*) , intent(in) :: cdfam ! family of observation
917 integer , intent(in) :: destObsColumn
918
919 ! Locals:
920 integer :: bodyIndex, headerIndex, levelIndex, bufrCode
921 real(8) :: observedDoppler, simulatedDoppler
922 real(8) :: levelAltLow, levelAltHigh
923 real(8) :: radarAltitude, beamAzimuth, beamElevation, obsAltitude
924 real(8) :: uuLow, uuHigh, vvLow, vvHigh, uuInterpolated, vvInterpolated
925 real(8) :: interpWeight
926
927 call obs_set_current_header_list(obsSpaceData, cdfam)
928 if (.not.beSilent) write(*,*) 'Entering subroutine oop_raDvel_nl, family: ', trim(cdfam)
929
930
931 !
932 ! Loop over all header indices of the 'RA' family with schema 'radvel':
933 !
934 HEADER: do
935
936 headerIndex = obs_getHeaderIndex(obsSpaceData)
937 if (headerIndex < 0) exit HEADER
938
939 radarAltitude = obs_headElem_r(obsSpaceData, OBS_ALT, headerIndex)
940 beamAzimuth = obs_headElem_r(obsSpaceData, OBS_RZAM, headerIndex)
941 beamElevation = obs_headElem_r(obsSpaceData, OBS_RELE, headerIndex)
942 call obs_set_current_body_list(obsSpaceData, headerIndex)
943 !
944 ! Loop over all body indices of the 'RA' family with schema 'radvel':
945 !
946 BODY: do
947 bodyIndex = obs_getBodyIndex(obsSpaceData)
948 if (bodyIndex < 0) exit BODY
949 ! Check that this observation has the expected bufr element ID
950 bufrCode = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex)
951 if (bufrCode /= bufr_radvel) cycle BODY
952 ! only process observations flagged to be assimilated
953 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated) cycle BODY
954
955 obsAltitude = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
956
957 ! Levels that bracket the observation from OBS_LYR
958 ! note to self: like in GEM, level=1 is the highest level
959 levelIndex = obs_bodyElem_i(obsSpaceData, OBS_LYR, bodyIndex)
960
961 levelAltHigh = col_getHeight(columnTrlOnTrlLev, levelIndex, headerIndex,'MM')
962 levelAltLow = col_getHeight(columnTrlOnTrlLev, levelIndex+1, headerIndex,'MM')
963
964 ! Vertical interpolation of model wind at observation height
965 interpWeight = (obsAltitude - levelAltLow)/(levelAltHigh - levelAltLow)
966 uuHigh = col_getElem(columnTrlOnTrlLev, levelIndex, headerIndex, 'UU')
967 uuLow = col_getElem(columnTrlOnTrlLev, levelIndex+1, headerIndex, 'UU')
968 vvHigh = col_getElem(columnTrlOnTrlLev, levelIndex, headerIndex, 'VV')
969 vvLow = col_getElem(columnTrlOnTrlLev, levelIndex+1, headerIndex, 'VV')
970 uuInterpolated = uuLow + interpWeight*(uuHigh - uuLow)
971 vvInterpolated = vvLow + interpWeight*(vvHigh - vvLow)
972
973 ! Doppler velocity is the projection of wind along direction of radar beam
974 ! Positive values indicates velocities "away" from the radar
975 simulatedDoppler = uuInterpolated*sin(beamAzimuth) + vvInterpolated*cos(beamAzimuth)
976
977 observedDoppler = obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex)
978
979 call obs_bodySet_r(obsSpaceData, destObsColumn, bodyIndex, observedDoppler-simulatedDoppler)
980
981 end do BODY
982 end do HEADER
983 if (.not. beSilent) write(*,*) 'Ending subroutine oop_raDvel_nl, family: ', trim(cdfam)
984
985 end subroutine oop_raDvel_nl
986
987 !--------------------------------------------------------------------------
988 ! oop_gpsro_nl
989 !--------------------------------------------------------------------------
990 subroutine oop_gpsro_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, destObsColumn)
991 !:Purpose: Computation of Jo and the residuals to the GPSRO observations
992 !
993 !:Note: gps_struct1sw_v2 allows calculation of partial derivatives of refractivity
994 ! in gps_diff object w.r.t TT/HU/GZ/P0. The indirect dependency refractivity
995 ! to TT/HU/P0 through GZ is now attributed to direct dependency of refractivity on GZ.
996 implicit none
997
998 ! Arguments:
999 type(struct_columnData) , intent(in) :: columnTrlOnTrlLev
1000 type(struct_obs) , intent(inout) :: obsSpaceData
1001 logical , intent(in) :: beSilent
1002 integer , intent(in) :: destObsColumn
1003
1004 ! Locals:
1005 type(struct_vco), pointer :: vco_hr
1006 real(8) :: pjob, pjo1
1007 real(8) :: zlat, lat
1008 real(8) :: zlon, lon
1009 real(8) :: zazm, azm
1010 integer :: isat, iclf
1011 real(8) :: rad, geo
1012 real(8), allocatable :: zpp(:)
1013 real(8), allocatable :: ztt(:)
1014 real(8), allocatable :: zhu(:)
1015 real(8), allocatable :: zHeight(:)
1016 real(8), allocatable :: zuu(:)
1017 real(8), allocatable :: zvv(:)
1018 real(8) :: zp0, zmt
1019 real(8) :: hnh1, zobs, zmhx, zoer, zinc
1020 integer headerIndex, idatyp, bodyIndex
1021 integer jl, ngpslev, nwndlev
1022 logical assim, firstheader, ldsc
1023 integer :: nh, nh1, iProfile, varNum
1024 type(gps_profile) :: prf
1025 real(8) , allocatable :: h (:),azmv(:)
1026 type(gps_diff), allocatable :: rstv(:)
1027 integer :: Vcode
1028
1029 if (.not.beSilent) write(*,*)'ENTER oop_gpsro_nl'
1030
1031 vco_hr => col_getVco(columnTrlOnTrlLev)
1032 vcode = vco_hr%vcode
1033
1034 !
1035 ! Initializations
1036 !
1037 ngpslev=col_getNumLev(columnTrlOnTrlLev,'TH')
1038 nwndlev=col_getNumLev(columnTrlOnTrlLev,'MM')
1039 allocate(zpp(ngpslev))
1040 allocate(ztt(ngpslev))
1041 allocate(zhu(ngpslev))
1042 allocate(zHeight(ngpslev))
1043 allocate(zuu(ngpslev))
1044 allocate(zvv(ngpslev))
1045
1046 allocate( h (gps_ro_maxprfsize) )
1047 allocate( azmv (gps_ro_maxprfsize) )
1048 allocate( rstv (gps_ro_maxprfsize) )
1049
1050 !
1051 ! Loop over all header indices of the 'RO' family:
1052 !
1053 call obs_set_current_header_list(obsSpaceData,'RO')
1054 firstheader = .true.
1055
1056 HEADER: do
1057 headerIndex = obs_getHeaderIndex(obsSpaceData)
1058 if (headerIndex < 0) exit HEADER
1059 !
1060 ! Process only refractivity data (codtyp 169)
1061 !
1062 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
1063 if ( idatyp /= 169 ) cycle HEADER
1064 !
1065 ! Scan for requested data values of the profile, and count them
1066 !
1067 assim = .false.
1068 nh = 0
1069 call obs_set_current_body_list(obsSpaceData, headerIndex)
1070 BODY: do
1071 bodyIndex = obs_getBodyIndex(obsSpaceData)
1072 if (bodyIndex < 0) exit BODY
1073 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
1074 assim = .true.
1075 nh = nh + 1
1076 end if
1077 end do BODY
1078 !
1079 ! If no assimilations are requested, skip to next header
1080 !
1081 if (.not.assim) cycle HEADER
1082 iProfile = gps_iprofile_from_index(headerIndex)
1083 varNum = gps_vRO_IndexPrf(iProfile, 2)
1084 !
1085 ! Basic geometric variables of the profile:
1086 !
1087 isat = obs_headElem_i(obsSpaceData,OBS_SAT,headerIndex)
1088 iclf = obs_headElem_i(obsSpaceData,OBS_ROQF,headerIndex)
1089 rad = obs_headElem_r(obsSpaceData,OBS_TRAD,headerIndex)
1090 geo = obs_headElem_r(obsSpaceData,OBS_GEOI,headerIndex)
1091 azm = obs_headElem_r(obsSpaceData,OBS_AZA,headerIndex)
1092 zmt = col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF')
1093 !
1094 ! Profile at the observation location:
1095 !
1096 zlat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
1097 zlon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
1098 lat = zlat * MPC_DEGREES_PER_RADIAN_R8
1099 lon = zlon * MPC_DEGREES_PER_RADIAN_R8
1100 zazm = azm / MPC_DEGREES_PER_RADIAN_R8
1101 zp0 = col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0')
1102 do jl = 1, ngpslev
1103 !
1104 ! Profile x
1105 !
1106 zpp(jl) = col_getPressure(columnTrlOnTrlLev,jl,headerIndex,'TH')
1107 ztt(jl) = col_getElem(columnTrlOnTrlLev,jl,headerIndex,'TT') - MPC_K_C_DEGREE_OFFSET_R8
1108 zhu(jl) = col_getElem(columnTrlOnTrlLev,jl,headerIndex,'HU')
1109 zHeight(jl) = col_getHeight(columnTrlOnTrlLev,jl,headerIndex,'TH')
1110 end do
1111
1112 if (Vcode == 5002) then
1113 ! case with top thermo level above top momentum level (Vcode=5002)
1114 do jl = 1, nwndlev
1115 zuu(jl) = col_getElem(columnTrlOnTrlLev,jl ,headerIndex,'UU')
1116 zvv(jl) = col_getElem(columnTrlOnTrlLev,jl ,headerIndex,'VV')
1117 end do
1118 elseif (Vcode == 5005 .or. Vcode == 5100 .or. Vcode == 21001) then
1119 ! case without top thermo above top momentum level or unstaggered (Vcode=5001/4/5)
1120 do jl = 1, nwndlev-1
1121 zuu(jl) = col_getElem(columnTrlOnTrlLev,jl+1,headerIndex,'UU')
1122 zvv(jl) = col_getElem(columnTrlOnTrlLev,jl+1,headerIndex,'VV')
1123 end do
1124 zuu(nwndlev) = zuu(nwndlev-1)
1125 zvv(nwndlev) = zuu(nwndlev-1)
1126 else
1127 ! case with Vcode /= 5002 and Vcode /= 5005 and Vcode /= 5100
1128 call utl_abort('oop_gpsro_nl: invalid vertical coord!')
1129 end if
1130 zuu(ngpslev) = zuu(nwndlev)
1131 zvv(ngpslev) = zuu(nwndlev)
1132 !
1133 ! GPS profile structure:
1134 !
1135 call gps_struct1sw_v2(ngpslev,zLat,zLon,zAzm,zMT,Rad,geo,zP0,zPP,zTT,zHU,zUU,zVV,zHeight,prf)
1136 ldsc=.not.btest(iclf,16-3)
1137 !
1138 ! Prepare the vector of all the observations:
1139 !
1140 nh1 = 0
1141 !
1142 ! Loop over all body indices for this headerIndex:
1143 ! (start at the beginning of the list)
1144 !
1145 call obs_set_current_body_list(obsSpaceData, headerIndex)
1146 BODY_2: do
1147 bodyIndex = obs_getBodyIndex(obsSpaceData)
1148 if (bodyIndex < 0) exit BODY_2
1149 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
1150 nh1 = nh1 + 1
1151 h(nh1) = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1152 azmv(nh1)= zazm
1153 end if
1154 end do BODY_2
1155 !
1156 ! Apply the observation operator:
1157 !
1158 ! varNum = bufr_nebd (15037) or varNum = bufr_nerf (15036) for GPS-RO
1159 iProfile = gps_iprofile_from_index(headerIndex)
1160 if (varNum == bufr_nebd) then
1161 call gps_bndopv1(h, azmv, nh, prf, rstv)
1162 else
1163 call gps_refopv (h, nh, prf, rstv)
1164 end if
1165 !
1166 ! Perform the (H(x)-Y)/S operation:
1167 !
1168 nh1 = 0
1169 pjob = 0.d0
1170 !
1171 ! Loop over all body indices for this headerIndex:
1172 ! (start at the beginning of the list)
1173 !
1174 call obs_set_current_body_list(obsSpaceData, headerIndex)
1175 BODY_3: do
1176 bodyIndex = obs_getBodyIndex(obsSpaceData)
1177 if (bodyIndex < 0) exit BODY_3
1178 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
1179 nh1 = nh1 + 1
1180 !
1181 ! Altitude:
1182 !
1183 hnh1= obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1184 if (varNum == bufr_nebd) hnh1=hnh1-rad
1185 !
1186 ! Observation operator H(x)
1187 !
1188 zmhx = rstv(nh1)%var
1189 !
1190 ! Observation value Y
1191 !
1192 zobs = obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
1193 !
1194 ! Observation error S
1195 !
1196 zoer = obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)
1197 !
1198 ! Normalized increment
1199 !
1200 zinc = (zmhx - zobs) / zoer
1201 !
1202 ! Datum contribution to Jo:
1203 !
1204 pjo1 = 0.5d0 * zinc * zinc
1205 !
1206 ! Per profile (PJOB) cumulatives:
1207 !
1208 pjob= pjob + pjo1
1209 !
1210 if (firstheader .and. .not.beSilent) then
1211 write(*, &
1212 '(A9,i10,3f7.2,f11.1,4f12.6,15f12.4)') 'DOBSGPSRO', &
1213 headerIndex,lat,lon,azm,hnh1,zobs,zoer, &
1214 zmhx,zinc,pjob,prf%gst(ngpslev)%var
1215 end if
1216 call obs_bodySet_r(obsSpaceData,destObsColumn,bodyIndex, zobs - zmhx)
1217 end if
1218 end do BODY_3
1219
1220 if ( .not.beSilent ) write(*,'(A9,i10,2f7.2,f18.10,f12.4,2I6)') &
1221 'GPSRO_JO',headerIndex,lat,lon,pjob,zmt,isat,ldsc
1222 firstheader = .false.
1223 end do HEADER
1224
1225 deallocate( rstv )
1226 deallocate( azmv )
1227 deallocate( h )
1228
1229 deallocate(zvv)
1230 deallocate(zuu)
1231 deallocate(zhu)
1232 deallocate(zHeight)
1233 deallocate(ztt)
1234 deallocate(zpp)
1235
1236 if (.not.beSilent) write(*,*)'EXIT oop_gpsro_nl'
1237
1238 end subroutine oop_gpsro_nl
1239
1240 !--------------------------------------------------------------------------
1241 ! oop_gpsgb_nl
1242 !--------------------------------------------------------------------------
1243 subroutine oop_gpsgb_nl( columnTrlOnTrlLev, obsSpaceData, beSilent, &
1244 destObsColumn, analysisMode_opt )
1245 !
1246 !:Purpose: Computation of the residuals to the GB-GPS ZTD observations
1247 !
1248 implicit none
1249
1250 ! Arguments:
1251 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
1252 type(struct_obs) , intent(inout) :: obsSpaceData
1253 logical , intent(in) :: beSilent
1254 integer , intent(in) :: destObsColumn
1255 logical, optional, intent(in) :: analysisMode_opt
1256
1257 ! Locals:
1258 real(8), allocatable :: zpp (:)
1259 real(8), allocatable :: ztt (:)
1260 real(8), allocatable :: zhu (:)
1261 real(8), allocatable :: zHeight (:)
1262 real(8) :: zlat, lat, zlon, lon
1263 real(8) :: zp0, zmt, zdzmin
1264 real(8) :: zobs, zoer, zinc, zhx, zlev
1265 real(8) :: zdz, zpsobs, zpsmod, zpwmod, zpomp, zpomps
1266 real(8) :: ztdomp(gps_gb_maxdata)
1267 real(8) :: bias, std
1268 integer :: headerIndex, bodyIndex, ioneobs, idatyp, bufrCode, index_ztd, iztd
1269 integer :: jl, nlev_T, nobs2p
1270 integer :: icount1, icount2, icount3, icount, icountp
1271 logical :: assim, llrej, analysisMode, lfsl
1272 character(len=12) :: cstnid
1273 type(gps_profilezd) :: prf
1274 type(gps_diff) :: ztdopv
1275 !
1276 ! PW lower limit (mm) and Ps O-P upper limit (Pa) for ZTD assimilation
1277 ! Note: 1 mb = 100 Pa --> 2.2 mm ZTD
1278 !
1279 real(8) :: zpwmin, zpompmx
1280 data zpwmin / 2.0d0 /
1281 data zpompmx / 200.0d0 /
1282 !
1283 ! Criteria to select single observation (1-OBS mode)
1284 !
1285 ! Minimum value for ZTD O-P (m)
1286 real(8) :: xompmin
1287 data xompmin / 0.015d0 /
1288 ! Minimum value for background (trial) PW (mm)
1289 real(8) :: xpwmin
1290 data xpwmin / 20.0d0 /
1291 ! Maximum height difference between observation and background surface (m)
1292 real(8) :: xdzmax
1293 data xdzmax / 400.0d0 /
1294
1295 if (.not.beSilent) write(*,*)'ENTER oop_gpsgb_nl'
1296
1297 if (present(analysisMode_opt)) then
1298 analysisMode = analysisMode_opt
1299 else
1300 analysisMode = .true.
1301 end if
1302
1303 zpomps = 0.0d0
1304
1305 ! Ensure Jacobian-related arrays are not allocated to force them to be recalculated in oop_H
1306 if (allocated(oop_vZTD_Jacobian)) deallocate(oop_vZTD_Jacobian)
1307
1308 zdzmin = gps_gb_dzmin
1309 nobs2p = 50
1310
1311 nlev_T = col_getNumLev(columnTrlOnTrlLev,'TH')
1312 if (gps_gb_ltestop .and. .not.beSilent) write(*,*) ' col_getNumLev[columnTrlOnTrlLev,TH] = ',nlev_T
1313
1314 !
1315 ! Initializations
1316 !
1317 allocate(ztt(nlev_T))
1318 allocate(zhu(nlev_T))
1319 allocate(zHeight(nlev_T))
1320 allocate(zpp(nlev_T))
1321
1322 if ( .not.beSilent ) then
1323 write(*, *) ' '
1324 write(*, *) ' '
1325 write(*,'(A11,A9,3A8,A9,4A8,2A9,A7,A10,A11)') &
1326 'OOP_GPSGB_NL','CSTNID','ZLAT','ZLON','ZLEV','ZDZ','ZOBS','ZOER','ZHX','O-P', &
1327 'ZPOMPS','ZPOMP','ZPWMOD','ZINC2'
1328 end if
1329
1330 icount = 0
1331 icount1 = 0
1332 icount2 = 0
1333 icount3 = 0
1334 icountp = 0
1335 ioneobs = -1
1336
1337 ! loop over all header indices of the 'GP' family (all obs locations/times)
1338 call obs_set_current_header_list(obsSpaceData,'GP')
1339
1340 HEADER: do
1341 headerIndex = obs_getHeaderIndex(obsSpaceData)
1342 if (headerIndex < 0) exit HEADER
1343
1344 ! Process only GP data (codtyp 189)
1345 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
1346 if ( idatyp /= 189 ) cycle HEADER
1347
1348 assim = .false.
1349 zpsobs = -100.0d0
1350 lfsl = .false.
1351 cstnid = obs_elem_c(obsSpaceData,'STID',headerIndex)
1352 if (index(cstnid,'FSL_') > 0 .or. index(cstnid,'-NOAA') > 0) lfsl = .true.
1353
1354 ! Scan for requested ZTD assimilation.
1355 ! Get GPS antenna height ZLEV and Ps(ZLEV) (ZPSOBS)
1356 !
1357 ! loop over all body indices for this headerIndex (observations at location/time)
1358 call obs_set_current_body_list(obsSpaceData, headerIndex)
1359 BODY: do
1360 bodyIndex = obs_getBodyIndex(obsSpaceData)
1361 if (bodyIndex < 0) exit BODY
1362 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
1363 if ( (bufrCode == bufr_nezd) .and. &
1364 (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) ) then
1365 zlev = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1366 assim = .true.
1367 ! Index in body of ZTD datum (assume at most 1 per header)
1368 index_ztd = bodyIndex
1369 icount = icount + 1
1370 end if
1371 if ( bufrCode == bufr_neps ) then
1372 if ( (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) .or. gps_gb_llblmet ) then
1373 zpsobs = obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
1374 zpomps = obs_bodyElem_r(obsSpaceData,destObsColumn,bodyIndex)
1375 end if
1376 end if
1377 end do BODY
1378
1379 ! If no ZTD assimilation requested, jump to next header
1380 if (.not.assim) cycle HEADER
1381
1382 ! Profile at the observation location:
1383 lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
1384 lon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
1385 zlat = lat * MPC_DEGREES_PER_RADIAN_R8
1386 zlon = lon * MPC_DEGREES_PER_RADIAN_R8
1387 zmt = col_getHeight(columnTrlOnTrlLev,0,headerIndex,'SF')
1388 zp0 = col_getElem(columnTrlOnTrlLev,1,headerIndex,'P0')
1389 do jl = 1, nlev_T
1390 zpp(jl) = col_getPressure(columnTrlOnTrlLev,jl,headerIndex,'TH')
1391 ztt(jl) = col_getElem(columnTrlOnTrlLev,jl,headerIndex,'TT')-MPC_K_C_DEGREE_OFFSET_R8
1392 zhu(jl) = col_getElem(columnTrlOnTrlLev,jl,headerIndex,'HU')
1393 zHeight(jl) = col_getHeight(columnTrlOnTrlLev,jl,headerIndex,'TH')
1394 end do
1395 zdz = zlev - zmt
1396
1397 ! Fill GPS ZTD profile structure (PRF):
1398 call gps_structztd_v2(nlev_T,lat,lon,zmt,zp0,zpp,ztt,zhu,zHeight,gps_gb_lbevis,gps_gb_irefopt,prf)
1399
1400 ! Apply the GPS ZTD observation operator
1401 ! --> output is model ZTD (type gps_diff) and P at obs height ZLEV
1402 call gps_ztdopv(zlev,prf,gps_gb_lbevis,zdzmin,ztdopv,zpsmod,gps_gb_iztdop)
1403
1404 ! Get model profile PW
1405 call gps_pw(prf,zpwmod)
1406 ! ZTD (m)
1407 zhx = ztdopv%var
1408
1409 ! If analysis mode, reject ZTD data for any of the following conditions:
1410 ! (1) the trial PW is too low (extremely dry)
1411 ! and if gps_gb_LASSMET=true and for NOAA/FSL sites only:
1412 ! (2) Ps observation is missing or out of normal range
1413 ! (3) the ABS(Ps(obs)-Ps(mod)) difference is too large
1414 llrej = .false.
1415 zpomp = -9999.0D0
1416 if ( analysisMode ) then
1417 llrej = ( zpwmod < zpwmin )
1418 if ( gps_gb_lassmet .and. lfsl ) then
1419 if ( .not. llrej ) then
1420 if ( zpsobs > 40000.0d0 .and. zpsobs <= 110000.0d0 ) then
1421 zpomp = zpsobs - zpsmod
1422 llrej = ( abs(zpomp) > zpompmx )
1423 if ( llrej ) icount3 = icount3 + 1
1424 else
1425 llrej = .true.
1426 icount2 = icount2 + 1
1427 end if
1428 else
1429 icount1 = icount1 + 1
1430 end if
1431 end if
1432 end if
1433
1434 if ( llrej ) then
1435 call obs_bodySet_i(obsSpaceData,OBS_ASS,index_ztd, obs_notAssimilated)
1436 if ( .not. gps_gb_lassmet ) icount1 = icount1 + 1
1437 end if
1438
1439 ! Perform the (H(x)-Y)/SDERR operation
1440 !
1441 ! loop over all body indices for this headerIndex
1442 call obs_set_current_body_list(obsSpaceData, headerIndex)
1443 BODY_2: do
1444 bodyIndex = obs_getBodyIndex(obsSpaceData)
1445 if (bodyIndex < 0) exit BODY_2
1446 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
1447 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated .and. &
1448 bufrCode == bufr_nezd ) then
1449 icountp = icountp + 1
1450 !
1451 ! Observation value Y
1452 !
1453 zobs = obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
1454 !
1455 ! Observation error SDERR
1456 !
1457 zoer = obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)
1458 if ( zoer <= 0.0d0 ) then
1459 write(*,*) ' Problem with ZTD observation error!'
1460 write(*,*) ' Station =',cstnid
1461 write(*,*) ' Error =', zoer
1462 call utl_abort('OOP_GPSGB_NL: ABORT! BAD ZTD OBSERR')
1463 end if
1464
1465 ! Observation height (m)
1466 !
1467 zlev = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1468 !
1469 ! Normalized increment ZINC
1470 !
1471 ztdomp(icountp) = zobs - zhx
1472 zinc = (zhx - zobs) / zoer
1473 call obs_bodySet_r(obsSpaceData,destObsColumn,bodyIndex, zobs - zhx)
1474
1475 !
1476 ! Apply data selection criteria for 1-OBS Mode
1477 !
1478 if ( gps_gb_l1obs .and. ioneobs == -1 ) then
1479 if ( (zobs-zhx) > xompmin .and. zpwmod > xpwmin .and. &
1480 abs(zdz) < xdzmax ) then
1481 ioneobs = headerIndex
1482 write(*,*) 'SINGLE OBS SITE = ',cstnid
1483 end if
1484 end if
1485 !
1486 ! Print data for first NOBS2P observations
1487 !
1488 if ( icountp <= nobs2p .and. .not.beSilent ) then
1489 write(*, &
1490 '(A14,A9,3(1x,f7.2),1x,f8.2,4(1x,f8.5),2(1x,f8.4),2x,f5.2,1x,f9.2,1x,f10.5)') &
1491 'OOP_GPSGB_NL: ',cstnid,zlat,zlon,zlev,zdz,zobs,zoer/gps_gb_yzderrwgt,zhx,-zinc*zoer, &
1492 zpomps/100.d0,zpomp/100.d0,zpwmod,zinc/zoer
1493 end if
1494
1495 end if
1496
1497 end do BODY_2
1498
1499 end do HEADER
1500
1501 deallocate(ztt)
1502 deallocate(zhu)
1503 deallocate(zHeight)
1504 deallocate(zpp)
1505
1506 if ( .not. beSilent ) then
1507 write(*,*) ' '
1508 write(*,*) 'NUMBER OF GPS ZTD DATA FLAGGED FOR ASSIMILATION = ', icountp
1509 if ( icountp > 0 ) then
1510 bias = sum(ztdomp(1:icountp))/real(icountp,8)
1511 std = 0.d0
1512 do jl = 1, icountp
1513 std = std + (ztdomp(jl)-bias)**2
1514 end do
1515 write(*, *) ' MEAN O-P (BIAS) [mm] = ', bias*1000.d0
1516 if (icountp > 1) then
1517 std = sqrt(std/(real(icountp,8)-1.d0))
1518 write(*, *) ' STD O-P [mm] = ', std*1000.d0
1519 else
1520 write(*, *) ' STD O-P Uncomputable since number of GPS ZTD observations is 1'
1521 end if
1522 write(*, *) ' '
1523 end if
1524 end if
1525
1526 if ( gps_gb_l1obs .and. analysisMode ) then
1527 ! Set assim flag to 0 for all observations except for selected record (site/time)
1528 if ( ioneobs /= -1 ) then
1529 call obs_set_current_header_list(obsSpaceData,'GP')
1530 icountp = 1
1531 HEADER_1: do
1532 headerIndex = obs_getHeaderIndex(obsSpaceData)
1533 if (headerIndex < 0) exit HEADER_1
1534 if (headerIndex /= ioneobs ) then
1535 call obs_set_current_body_list(obsSpaceData, headerIndex)
1536 BODY_1: do
1537 bodyIndex = obs_getBodyIndex(obsSpaceData)
1538 if (bodyIndex < 0) exit BODY_1
1539 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex, obs_notAssimilated)
1540 end do BODY_1
1541 end if
1542 end do HEADER_1
1543 else
1544 call utl_abort('ERROR: FAILED TO SELECT SINGLE OBSERVATION!')
1545 end if
1546 end if
1547
1548 gps_gb_numztd = icountp
1549
1550 if ( analysisMode .and. icount > 0 .and. .not.gps_gb_l1obs .and. .not.beSilent ) then
1551 write(*,*) ' '
1552 write(*,*) '-----------------------------------------'
1553 write(*,*) ' SUMMARY OF ZTD REJECTIONS IN OOP_GPSGB_NL '
1554 write(*,*) '-----------------------------------------'
1555 write(*,*) ' TOTAL NUMBER OF ZTD DATA ORIGINALLY FLAGGED FOR ASSMILATION = ', icount
1556 write(*,*) ' NUMBER OF ZTD DATA REJECTED DUE TO LOW TRIAL PW = ', icount1
1557 write(*,*) ' NUMBER OF ZTD DATA REJECETD DUE TO NO PS OBS = ', icount2
1558 write(*,*) ' NUMBER OF ZTD DATA REJECETD DUE TO LARGE PS O-P = ', icount3
1559 write(*,*) ' TOTAL NUMBER OF REJECTED ZTD DATA = ', icount1+icount2+icount3
1560 write(*,*) ' PERCENT REJECTED = ', &
1561 (real(icount1+icount2+icount3,8) / real(icount,8))*100.0d0
1562 write(*, *) ' TOTAL NUMBER OF ASSIMILATED ZTD DATA = ', icountp
1563 write(*,*) ' '
1564 end if
1565
1566 if ( icount > 0 .and. gps_gb_numztd > 0) then
1567
1568 if ( analysisMode ) then
1569 if ( .not.beSilent ) write(*,*) ' Number of GPS ZTD data to be assimilated (gps_gb_numZTD) = ', gps_gb_numztd
1570 else
1571 if ( .not.beSilent ) write(*,*) ' Number of GPS ZTD data for background check (gps_gb_numZTD) = ', gps_gb_numztd
1572 end if
1573
1574 if ( .not.beSilent ) write(*,*) ' Allocating and setting gps_ZTD_Index(gps_gb_numZTD)...'
1575 if (allocated(gps_ZTD_Index)) deallocate(gps_ZTD_Index)
1576 allocate(gps_ZTD_Index(gps_gb_numztd))
1577 iztd = 0
1578 call obs_set_current_header_list(obsSpaceData,'GP')
1579 HEADER_2: do
1580 headerIndex = obs_getHeaderIndex(obsSpaceData)
1581 if (headerIndex < 0) exit HEADER_2
1582 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
1583 if ( idatyp == 189 ) then
1584 call obs_set_current_body_list(obsSpaceData, headerIndex)
1585 BODY_3: do
1586 bodyIndex = obs_getBodyIndex(obsSpaceData)
1587 if (bodyIndex < 0) exit BODY_3
1588 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
1589 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated .and. &
1590 bufrCode == bufr_nezd ) then
1591 iztd = iztd + 1
1592 gps_ZTD_Index(iztd) = headerIndex
1593 end if
1594 end do BODY_3
1595 end if
1596 end do HEADER_2
1597
1598 if ( iztd /= gps_gb_numztd ) then
1599 call utl_abort('ERROR: gps_ZTD_Index init: iztd /= gps_gb_numZTD!')
1600 end if
1601
1602 end if
1603
1604 if (.not.beSilent) write(*,*)'EXIT oop_gpsgb_nl'
1605
1606 end subroutine oop_gpsgb_nl
1607
1608 !--------------------------------------------------------------------------
1609 ! oop_tovs_nl
1610 !--------------------------------------------------------------------------
1611 subroutine oop_tovs_nl( columnTrl, obsSpaceData, datestamp, beSilent, &
1612 bgckMode_opt, option_opt, sourceObs_opt, destObs_opt )
1613 !:Purpose: Computation of the residuals to the tovs observations
1614 ! option_opt: defines input state:
1615 ! - 'HR': High Resolution background state,
1616 ! - 'LR': Low Resolution background state, (CURRENTLY NOT SUPPORTED)
1617 ! - 'MO': Model state. (CURRENTLY NOT SUPPORTED)
1618 implicit none
1619
1620 ! Arguments:
1621 type(struct_columnData), intent(in) :: columnTrl
1622 type(struct_obs), intent(inout) :: obsSpaceData
1623 integer, intent(in) :: datestamp
1624 logical, intent(in) :: beSilent
1625 logical, optional, intent(in) :: bgckMode_opt
1626 character(len=*), optional, intent(in) :: option_opt ! only valid value is HR
1627 integer, optional, intent(in) :: sourceObs_opt ! usually set to OBS_VAR
1628 integer, optional, intent(in) :: destObs_opt ! usually set to OBS_OMP
1629
1630 ! Locals:
1631 integer :: jdata, sourceObs, destObs
1632 logical :: llprint,bgckMode
1633 character(len=2) :: option
1634 integer :: channelIndex, tovsIndex
1635 real(pre_obsReal) :: zdtb, obsPRM
1636 integer :: idatyp, channelNumber
1637 integer :: headerIndex, bodyIndex
1638
1639 if (.not.obs_famExist(obsSpaceData,'TO', localMPI_opt = .true. )) return
1640
1641 ! 0. set default values if bgckMode, option and source/dest columns not specified
1642 !
1643
1644 if (.not.beSilent) write(*,*) "Entering subroutine oop_tovs_nl"
1645
1646 if (present(bgckMode_opt)) then
1647 bgckMode = bgckMode_opt
1648 else
1649 bgckMode = .false.
1650 end if
1651
1652 if (present(option_opt)) then
1653 option = option_opt(1:2)
1654 else
1655 option = 'HR'
1656 end if
1657 if ( option /= 'HR' ) call utl_abort('oop_tovs_nl: Invalid option for input state')
1658
1659 if (present(sourceObs_opt)) then
1660 sourceObs = sourceObs_opt
1661 else
1662 sourceObs = OBS_VAR
1663 end if
1664
1665 if (present(destObs_opt)) then
1666 destObs = destObs_opt
1667 else
1668 destObs = OBS_OMP
1669 end if
1670
1671 ! 1. Prepare atmospheric profiles for all tovs observation points for use in rttov
1672 ! . -----------------------------------------------------------------------------
1673 call tvs_fillProfiles(columnTrl,obsSpaceData,datestamp,"nl",beSilent)
1674
1675 if ( .not.beSilent ) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1676
1677 ! 2. Compute radiance
1678 ! . ----------------
1679 call tvs_rttov(obsSpaceData,bgckMode,beSilent)
1680 if ( .not.beSilent ) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1681
1682 ! 3. Compute the residuals
1683 ! . ----------------------------
1684 if ( option == 'HR' .or. option == 'LR' ) then
1685 do jdata=1,obs_numbody(obsSpaceData)
1686 call obs_bodySet_r(obsSpaceData,OBS_PRM,jdata, obs_bodyElem_r(obsSpaceData,sourceObs,jdata))
1687 end do
1688 end if
1689
1690 ! loop over all header indices of the 'TO' family
1691 call obs_set_current_header_list(obsSpaceData,'TO')
1692 HEADER: do
1693 headerIndex = obs_getHeaderIndex(obsSpaceData)
1694 if (headerIndex < 0) exit HEADER
1695
1696 ! Extract general information for this observation point
1697 ! ------------------------------------------------------
1698
1699 ! process only radiance data to be assimilated?
1700 idatyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
1701 if ( .not. tvs_isIdBurpTovs(idatyp) ) then
1702 write(*,*) 'oop_tovs_nl: warning unknown radiance codtyp present check NAMTOVSINST', idatyp
1703 cycle HEADER
1704 end if
1705 tovsIndex = tvs_tovsIndex(headerIndex)
1706 if ( tovsIndex == -1 ) cycle HEADER
1707
1708 ! Set the body list
1709 ! (& start at the beginning of the list)
1710 call obs_set_current_body_list(obsSpaceData, headerIndex)
1711 BODY: do
1712 bodyIndex = obs_getBodyIndex(obsSpaceData)
1713 if ( bodyIndex < 0 ) exit BODY
1714
1715 ! Only consider if flagged for assimilation
1716 if ( obs_bodyElem_i(obsSpaceData,obs_ASS,bodyIndex) /= obs_assimilated ) cycle BODY
1717
1718 call tvs_getChannelNumIndexFromPPP( obsSpaceData, headerIndex, bodyIndex, &
1719 channelNumber, channelIndex )
1720
1721 if ( channelIndex == 0 ) call utl_abort('oop_tovs_nl: error with channel number')
1722
1723 zdtb = obs_bodyElem_r(obsSpaceData,OBS_PRM,bodyIndex) - &
1724 tvs_radiance (tovsIndex) % bt(channelIndex)
1725 if ( tvs_debug ) then
1726 obsPRM = obs_bodyElem_r(obsSpaceData,OBS_PRM,bodyIndex)
1727 write(*,'(a,i4,2f8.2,f6.2)') ' channelNumber,sim,obs,diff= ', &
1728 channelNumber, tvs_radiance (tovsIndex) % bt(channelIndex), &
1729 obsPRM, -zdtb
1730 end if
1731 call obs_bodySet_r(obsSpaceData,destObs,bodyIndex, zdtb)
1732
1733 ! inflate OBS_OER for all-sky assimilation
1734 call oer_inflateErrAllsky(obsSpaceData, bodyIndex, destObs, beSilent_opt=.true.)
1735
1736 end do BODY
1737
1738 end do HEADER
1739
1740 if (option == 'HR') then
1741 llprint = .true.
1742 else
1743 llprint = .false.
1744 end if
1745 if ( beSilent ) llprint = .false.
1746 if ( llprint ) call tvs_printDetailledOmfStatistics(obsSpaceData)
1747
1748 end subroutine oop_tovs_nl
1749
1750 !--------------------------------------------------------------------------
1751 ! oop_chm_nl
1752 !--------------------------------------------------------------------------
1753 subroutine oop_chm_nl( columnTrlOnTrlLev, obsSpaceData, destObsColumn )
1754 !:Purpose: Computation of the residuals to the observations
1755 ! for all observations of the CH (chemical constituents) family.
1756 ! The array columnTrlOnTrlLev contains the input model array.
1757 ! Stores OmP in OBS_OMP in obsSpaceData.
1758 implicit none
1759
1760 ! Arguments:
1761 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
1762 type(struct_obs) , intent(inout) :: obsSpaceData
1763 integer , intent(in) :: destObsColumn
1764
1765 if (.not.obs_famExist(obsSpaceData,'CH', localMPI_opt = .true. )) return
1766
1767 if ( destObsColumn /= obs_omp ) then
1768 write(*,*) 'oop_chm_nl: WARNING: Storing results in an obs column other than OBS_OMP. Not fully implemented.'
1769 end if
1770
1771 call oopc_CHobsoperators( columnTrlOnTrlLev,obsSpaceData,'nl', & ! 'nl' for non-linear operator
1772 destObsColumn_opt=destObsColumn )
1773
1774 end subroutine oop_chm_nl
1775
1776 !--------------------------------------------------------------------------
1777 ! oop_Htl
1778 !--------------------------------------------------------------------------
1779 subroutine oop_Htl( columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData, &
1780 min_nsim, initializeLinearization_opt )
1781 !
1782 !:Purpose: Compute simulated observations from profiled model increments.
1783 ! It returns Hdx in OBS_WORK. Calls the several linear observation operators.
1784 implicit none
1785
1786 ! Arguments:
1787 type(struct_columnData), intent(in) :: columnAnlInc
1788 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev
1789 type(struct_obs) , intent(inout) :: obsSpaceData
1790 integer, intent(in) :: min_nsim
1791 logical, optional, intent(in) :: initializeLinearization_opt
1792
1793 ! Locals:
1794 type(struct_vco), pointer :: vco_anl
1795 logical, save :: initializeLinearization = .true.
1796
1797 if ( mmpi_myid == 0 ) then
1798 write(*,*) 'OOP_Htl - Linearized observation operators'
1799 end if
1800
1801 vco_anl => col_getVco(columnTrlOnAnlIncLev)
1802
1803 ! Re-linearlize if it is asked for
1804 if ( present(initializeLinearization_opt) ) then
1805 initializeLinearization = initializeLinearization_opt
1806 end if
1807
1808 if ( initializeLinearization ) then
1809 ! Find interpolation layer in model profiles (used by several operators)
1810 if ( col_getNumLev(columnTrlOnAnlIncLev,'MM') > 1 ) call oop_vobslyrs(columnTrlOnAnlIncLev, obsSpaceData, beSilent=.false.)
1811
1812 ! Initialize some operators needed by linearized H
1813 call subasic_obs(columnTrlOnAnlIncLev)
1814
1815 initializeLinearization = .false.
1816 end if
1817
1818 call oop_Hpp() ! fill in OBS_WORK : Hdx
1819
1820 call oop_Hsf() ! fill in OBS_WORK : Hdx
1821
1822 call oop_Hto() ! fill in OBS_WORK : Hdx
1823
1824 call oop_Hro( initializeLinearization_opt=initializeLinearization_opt )
1825
1826 if ( gps_gb_numZTD > 0 ) then
1827 call oop_Hgp( initializeLinearization_opt=initializeLinearization_opt )
1828 end if
1829
1830 call oop_Hchm() ! fill in OBS_WORK : Hdx
1831
1832 call oop_Hsst() ! fill in OBS_WORK : Hdx
1833
1834 call oop_Hice() ! fill in OBS_WORK : Hdx
1835
1836 call oop_Hhydro() ! fill in OBS_WORK : Hdx
1837
1838 call oop_HheightCoordObs() ! fill in OBS_WORK : Hdx
1839
1840 contains
1841
1842 !--------------------------------------------------------------------------
1843
1844 subroutine subasic_obs( columnTrlOnAnlIncLev )
1845 !:Purpose: Initialise background state dependant factors
1846 ! and vectors for use in TLM and adjoint of
1847 ! non-linear operator
1848 implicit none
1849
1850 ! Arguments:
1851 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev
1852
1853 ! Locals:
1854 type(struct_vco), pointer :: vco_anl
1855 integer :: jlev,columnIndex,nlev_T,vcode_anl,status
1856 real(8) :: zhu,one
1857
1858 if ( .not.col_varExist(columnTrlOnAnlIncLev,'TT') .or. .not.col_varExist(columnTrlOnAnlIncLev,'HU') ) return
1859
1860 write(*,*) 'subasic_obs: setting up linearized Tv operator'
1861
1862 vco_anl => col_getVco(columnTrlOnAnlIncLev)
1863 one=1.0D0
1864 nlev_T = col_getNumLev(columnTrlOnAnlIncLev,'TH')
1865 Vcode_anl = vco_anl%vCode
1866
1867 if ( Vcode_anl /= 5002 .and. Vcode_anl /= 5005 ) then
1868 call utl_abort('subasic_obs: invalid vertical coord!')
1869 end if
1870
1871 ! initialize virtual temperature operator
1872
1873 !$OMP PARALLEL DO PRIVATE(jlev,columnIndex,zhu)
1874 do jlev = 1, nlev_T
1875 do columnIndex=1,col_getNumCol(columnTrlOnAnlIncLev)
1876 zhu=col_getElem(columnTrlOnAnlIncLev,jlev,columnIndex,'HU')
1877 columnTrlOnAnlIncLev%oltv(1,jlev,columnIndex) = phf_fottva(zhu,one)
1878 columnTrlOnAnlIncLev%oltv(2,jlev,columnIndex) = phf_folnqva(zhu,col_getElem(columnTrlOnAnlIncLev, &
1879 jlev,columnIndex,'TT'),one)
1880 end do
1881 end do
1882 !$OMP END PARALLEL DO
1883
1884 end subroutine subasic_obs
1885
1886 !--------------------------------------------------------------------------
1887
1888 subroutine oop_Hpp()
1889 !: Purpose: Compute simulated Upper Air observations from profiled model
1890 ! increments.
1891 ! It returns Hdx in OBS_WORK
1892 ! Interpolate vertically the contents of commvo to
1893 ! the pressure levels of the observations.
1894 ! A linear interpolation in ln(p) is performed.
1895 implicit none
1896
1897 ! Locals:
1898 integer levIndexBot,levIndexTop
1899 integer headerIndex,INDEX_FAMILY,layerIndex
1900 integer J,bodyIndex,bufrCode,nlev_T
1901 real(8) ZDADPS
1902 real(8) ZWB,ZWT,ZLTV
1903 real(8) ZLEV,ZPT,ZPB
1904 real(8) delPT,delPB
1905 real(8) anlIncValueBot,anlIncValueTop,trlValueBot,trlValueTop
1906 integer, parameter :: numFamily=3
1907 character(len=2) :: list_family(numFamily)
1908 character(len=4) :: varLevel
1909
1910 list_family(1) = 'UA'
1911 list_family(2) = 'AI'
1912 list_family(3) = 'SW'
1913
1914 FAMILY: do index_family=1,numFamily
1915
1916 call obs_set_current_body_list(obsSpaceData,list_family(index_family))
1917 BODY: do
1918 bodyIndex = obs_getBodyIndex(obsSpaceData)
1919 if (bodyIndex < 0) exit BODY
1920
1921 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated .and. &
1922 obs_bodyElem_i(obsSpaceData,OBS_XTR,bodyIndex) == 0 .and. &
1923 obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) == 2 ) then
1924 headerIndex = obs_bodyElem_i(obsSpaceData,OBS_HIND,bodyIndex)
1925 ZLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1926 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
1927 varLevel = vnl_varLevelFromVarnum(bufrCode)
1928 layerIndex = obs_bodyElem_i(obsSpaceData,OBS_LYR,bodyIndex)
1929 levIndexTop = layerIndex + col_getOffsetFromVarno(columnTrlOnAnlIncLev,bufrCode)
1930 levIndexBot = levIndexTop+1
1931 ZPT = col_getPressure(columnTrlOnAnlIncLev,layerIndex ,headerIndex,varLevel)
1932 ZPB = col_getPressure(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,varLevel)
1933 delPT = col_getPressure(columnAnlInc,layerIndex ,headerIndex,varLevel)
1934 delPB = col_getPressure(columnAnlInc,layerIndex+1,headerIndex,varLevel)
1935 ZWB = LOG(ZLEV/ZPT)/LOG(ZPB/ZPT)
1936 ZWT = 1.0D0 - ZWB
1937
1938 ZDADPS = ( LOG(ZLEV/ZPB)*delPT/ZPT - &
1939 LOG(ZLEV/ZPT)*delPB/ZPB )/ &
1940 LOG(ZPB/ZPT)**2
1941
1942 if ( bufrCode == bufr_nees ) then
1943 anlIncValueBot=hutoes_tl(col_getElem(columnAnlInc,layerIndex+1,headerIndex,'HU'), &
1944 col_getElem(columnAnlInc,layerIndex+1,headerIndex,'TT'), &
1945 delPB, &
1946 col_getElem(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'HU'), &
1947 col_getPressure(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'TH'))
1948 anlIncValueTop=hutoes_tl(col_getElem(columnAnlInc,layerIndex ,headerIndex,'HU'), &
1949 col_getElem(columnAnlInc,layerIndex ,headerIndex,'TT'), &
1950 delPT, &
1951 col_getElem(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'HU'), &
1952 col_getPressure(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'TH'))
1953 trlValueBot=hutoes(col_getElem(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'HU'), &
1954 col_getElem(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'TT'), &
1955 col_getPressure(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'TH'))
1956 trlValueTop=hutoes(col_getElem(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'HU'), &
1957 col_getElem(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'TT'), &
1958 col_getPressure(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'TH'))
1959 else
1960 anlIncValueBot=col_getElem(columnAnlInc,levIndexBot,headerIndex)
1961 anlIncValueTop=col_getElem(columnAnlInc,levIndexTop,headerIndex)
1962 trlValueBot=col_getElem(columnTrlOnAnlIncLev,levIndexBot,headerIndex)
1963 trlValueTop=col_getElem(columnTrlOnAnlIncLev,levIndexTop,headerIndex)
1964 end if
1965 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex, &
1966 ZWB*anlIncValueBot + ZWT*anlIncValueTop+ &
1967 (trlValueBot - trlValueTop)* &
1968 ZDADPS)
1969
1970 end if
1971
1972 end do BODY
1973
1974 end do FAMILY
1975
1976 end subroutine oop_Hpp
1977
1978 !--------------------------------------------------------------------------
1979
1980 subroutine oop_Hsf()
1981 !:Purpose: Compute simulated surface observations from profiled model
1982 ! increments.
1983 ! It returns Hdx in OBS_WORK
1984 !
1985 implicit none
1986
1987 ! Locals:
1988 integer :: levIndexBot,levIndexTop
1989 integer :: headerIndex,layerIndex
1990 integer :: J, bodyIndex, bufrCode, INDEX_FAMILY, nlev, nLev_T
1991 real(8) :: coeffA, coeffB
1992 real(8) :: ZWB, ZWT, ZLTV, trlVirtTemp
1993 real(8) :: ZPT, ZPB, ZDELPS, ZDELTV, deltaT, obsHeight
1994 real(8) :: anlIncValue
1995 real(8) :: delP
1996 real(8) :: anlIncUwind, anlIncVwind, trlUwind, trlVwind, squareSum, trlWindSpeed, anlIncWindSpeed
1997 integer, parameter :: numFamily=5
1998 character(len=2) :: list_family(numFamily)
1999 character(len=4) :: varLevel
2000
2001 ! Temperature lapse rate for extrapolation of height below model surface
2002 coeffB = 1.0d0/(MPC_RGAS_DRY_AIR_R8*temperatureLapseRate/ec_rg)
2003 !
2004 !
2005 list_family(1) = 'UA'
2006 list_family(2) = 'SF'
2007 list_family(3) = 'SC'
2008 list_family(4) = 'GP'
2009 list_family(5) = 'RA'
2010
2011 FAMILY: do index_family=1,numFamily
2012
2013 call obs_set_current_body_list(obsSpaceData, list_family(index_family))
2014 BODY: do
2015 bodyIndex = obs_getBodyIndex(obsSpaceData)
2016 if (bodyIndex < 0) exit BODY
2017
2018 ! Process all data within the domain of the model
2019 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
2020 if ( bufrCode == bufr_nezd .or. bufrCode == bufr_radvel ) cycle BODY
2021 if ( (obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) == 1) &
2022 .and. (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) &
2023 .and. (bufrCode == bufr_nets .or. bufrCode == bufr_neps &
2024 .or. bufrCode == bufr_nepn .or. bufrCode == bufr_ness &
2025 .or. bufrCode == bufr_neus .or. bufrCode == bufr_nevs &
2026 .or. bufrCode == bufr_vis .or. bufrCode == bufr_logVis &
2027 .or. bufrCode == bufr_gust .or. bufrCode == bufr_nefs &
2028 .or. bufrCode == bufr_radarPrecip .or. bufrCode == bufr_logRadarPrecip &
2029 .or. obs_bodyElem_i(obsSpaceData,OBS_XTR,bodyIndex) == 0) ) then
2030
2031 varLevel = vnl_varLevelFromVarnum(bufrCode)
2032 nlev = col_getNumLev(columnAnlInc,varLevel)
2033 nlev_T = col_getNumLev(columnAnlInc,'TH')
2034 headerIndex = obs_bodyElem_i(obsSpaceData,OBS_HIND,bodyIndex)
2035 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
2036 layerIndex = obs_bodyElem_i(obsSpaceData,OBS_LYR,bodyIndex)
2037 obsHeight = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
2038 levIndexTop = nlev - 1 + col_getOffsetFromVarno(columnTrlOnAnlIncLev,bufrCode)
2039 levIndexBot = levIndexTop + 1
2040
2041 if (bufrCode == bufr_nets .OR. bufrCode == bufr_ness .OR. &
2042 bufrCode == bufr_neus .OR. bufrCode == bufr_nevs .OR. &
2043 bufrCode == bufr_vis .or. bufrCode == bufr_logVis .or. &
2044 bufrCode == bufr_gust .or. &
2045 bufrCode == bufr_radarPrecip .or. bufrCode == bufr_logRadarPrecip) THEN
2046 if (bufrCode == bufr_ness ) THEN
2047 delP = col_getPressure(columnAnlInc,nlev_T,headerIndex,'TH')
2048 anlIncValue = hutoes_tl(col_getElem(columnAnlInc,nlev_T,headerIndex,'HU'), &
2049 col_getElem(columnAnlInc,nlev_T,headerIndex,'TT'), &
2050 delP, &
2051 col_getElem(columnTrlOnAnlIncLev,nlev_T,headerIndex,'HU'), &
2052 col_getPressure(columnTrlOnAnlIncLev,nlev,headerIndex,varLevel))
2053 else
2054 anlIncValue=col_getElem(columnAnlInc,levIndexBot,headerIndex)
2055 end if
2056 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex,anlIncValue)
2057 else if (bufrCode == bufr_nefs) then
2058 anlIncUwind=col_getElem(columnAnlInc,col_getNumLev(columnAnlInc,varLevel),headerIndex,'UU')
2059 anlIncVwind=col_getElem(columnAnlInc,col_getNumLev(columnAnlInc,varLevel),headerIndex,'VV')
2060 trlUwind=col_getElem(columnTrlOnAnlIncLev,col_getNumLev(columnTrlOnAnlIncLev,varLevel),headerIndex,'UU')
2061 trlVwind=col_getElem(columnTrlOnAnlIncLev,col_getNumLev(columnTrlOnAnlIncLev,varLevel),headerIndex,'VV')
2062 squareSum=trlUwind**2+trlVwind**2
2063 if ( squareSum .gt. 1.d-10 ) then
2064 trlWindSpeed=sqrt(squareSum)
2065 anlIncWindSpeed=( trlUwind * anlIncUwind + trlVwind * anlIncVwind ) / trlWindSpeed
2066 else
2067 trlWindSpeed=0.0
2068 anlIncWindSpeed=0.0
2069 end if
2070 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex,anlIncWindSpeed)
2071 else if (bufrCode == bufr_neps .OR. bufrCode == bufr_nepn) THEN
2072 ZLTV = columnTrlOnAnlIncLev%OLTV(1,nlev_T,headerIndex)*col_getElem(columnAnlInc,nlev_T,headerIndex,'TT') &
2073 + columnTrlOnAnlIncLev%OLTV(2,nlev_T,headerIndex)*col_getElem(columnAnlInc,nlev_T,headerIndex,'HU')
2074 trlVirtTemp = columnTrlOnAnlIncLev%OLTV(1,nlev_T,headerIndex)*col_getElem(columnTrlOnAnlIncLev,nlev_T,headerIndex,'TT')
2075 deltaT= temperatureLapseRate*(obsHeight-col_getHeight(columnTrlOnAnlIncLev,0,headerIndex,'SF'))
2076 coeffA = ((trlVirtTemp-deltaT)/trlVirtTemp)
2077 ZDELPS= (col_getElem(columnAnlInc,1,headerIndex,'P0')*coeffA**coeffB)
2078 ZDELTV= ((col_getElem(columnTrlOnAnlIncLev,1,headerIndex,'P0')*coeffB*coeffA**(coeffB-1)) &
2079 *(deltaT/(trlVirtTemp*trlVirtTemp)*ZLTV))
2080 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex, ZDELPS+ZDELTV)
2081 else
2082
2083 call utl_abort('oop_Hsf: You have entered the twilight zone!')
2084
2085 end if
2086
2087 end if
2088
2089 end do BODY
2090
2091 end do FAMILY
2092
2093 end subroutine oop_Hsf
2094
2095 !--------------------------------------------------------------------------
2096
2097 subroutine oop_Hsst()
2098 !:Purpose: Compute simulated sea surface temperature observations
2099 ! from profiled model increments.
2100 ! It returns Hdx in OBS_WORK
2101 implicit none
2102
2103 ! Locals:
2104 integer :: headerIndex, bodyIndex, bufrCode
2105 real(8) :: anlIncValueBot
2106 character(len=4) :: varName
2107
2108 call obs_set_current_body_list( obsSpaceData, 'TM' )
2109
2110 !$OMP PARALLEL DO PRIVATE(bodyIndex, bufrCode, varName, headerIndex, anlIncValueBot)
2111 BODY: do bodyIndex = 1, obs_numBody( obsSpaceData )
2112 if ( obs_getFamily(obsSpaceData, bodyIndex_opt=bodyIndex) /= 'TM' ) cycle BODY
2113
2114 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
2115 if ( bufrCode /= bufr_sst ) cycle BODY
2116
2117 if ( col_varExist(columnAnlInc,'TM') ) then
2118 varName = 'TM'
2119 else
2120 varName = 'TG'
2121 end if
2122
2123 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated ) then
2124
2125 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
2126 anlIncValueBot = col_getElem(columnAnlInc, 1, headerIndex, varName_opt = varName )
2127 call obs_bodySet_r( obsSpaceData, OBS_WORK, bodyIndex, anlIncValueBot )
2128 end if
2129
2130 end do BODY
2131 !$OMP END PARALLEL DO
2132
2133 end subroutine oop_Hsst
2134
2135 !--------------------------------------------------------------------------
2136
2137 subroutine oop_Hhydro()
2138 !:Purpose: Compute simulated hydrological observations
2139 ! from profiled model increments.
2140 ! It returns Hdx in OBS_WORK
2141 implicit none
2142
2143 ! Locals:
2144 integer :: headerIndex, bodyIndex, bufrCode
2145 real(8) :: anlIncValueBot
2146 character(len=4) :: varName
2147
2148 call obs_set_current_body_list( obsSpaceData, 'HY' )
2149
2150 BODY: do
2151 bodyIndex = obs_getBodyIndex( obsSpaceData )
2152 if (bodyIndex < 0) exit BODY
2153
2154 ! Process all data within the domain of the model
2155 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
2156
2157 if ( bufrCode /= bufr_riverFlow ) cycle BODY
2158
2159 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated ) then
2160 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
2161 varName = vnl_varNameFromVarNum(bufrCode)
2162 anlIncValueBot = col_getElem(columnAnlInc, 1, headerIndex,varName_opt = varName)
2163 call obs_bodySet_r( obsSpaceData, OBS_WORK, bodyIndex, anlIncValueBot )
2164 end if
2165
2166 end do BODY
2167
2168 end subroutine oop_Hhydro
2169
2170 !--------------------------------------------------------------------------
2171
2172 subroutine oop_Hice()
2173 !:Purpose: Compute simulated sea ice concentration observations
2174 ! from profiled model increments.
2175 ! It returns Hdx in OBS_WORK
2176 implicit none
2177
2178 ! Locals:
2179 integer :: headerIndex, bodyIndex, bufrCode
2180 integer :: idate, imonth
2181 integer :: trackCellNum
2182 real(8) :: anlIncValue, scaling
2183 character(len=4) :: varName
2184
2185 call obs_set_current_body_list( obsSpaceData, 'GL' )
2186
2187 BODY: do
2188
2189 bodyIndex = obs_getBodyIndex( obsSpaceData )
2190 if (bodyIndex < 0) exit BODY
2191
2192 ! Process all data within the domain of the model
2193 scaling = oop_iceScaling(obsSpaceData, bodyIndex)
2194
2195 if (scaling == 0.0d0) cycle BODY
2196
2197 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated ) then
2198 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
2199 varName = vnl_varNameFromVarNum(bufrCode)
2200 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
2201 anlIncValue = scaling*col_getElem( columnAnlInc, 1, headerIndex, varName_opt = varName )
2202 call obs_bodySet_r( obsSpaceData, OBS_WORK, bodyIndex, anlIncValue )
2203 end if
2204
2205 end do BODY
2206
2207 end subroutine oop_Hice
2208
2209 !--------------------------------------------------------------------------
2210
2211 subroutine oop_Hto()
2212 !:Purpose: Compute simulated radiances observations from profiled model
2213 ! increments.
2214 ! It returns Hdx in OBS_WORK
2215 implicit none
2216
2217 ! Locals:
2218 integer :: datestamp
2219
2220 if (.not.obs_famExist(obsSpaceData,'TO', localMPI_opt = .true. )) return
2221
2222 ! 1. Prepare atmospheric profiles for all tovs observation points for use in rttov
2223 ! . -----------------------------------------------------------------------------
2224 !
2225 if (min_nsim == 1) then
2226 datestamp = tim_getDatestamp()
2227 call tvs_fillProfiles(columnTrlOnAnlIncLev, obsSpaceData, datestamp, "tlad", .false.)
2228 end if
2229
2230 ! 2. Compute radiance
2231 ! . ----------------
2232 !
2233 call tvslin_rttov_tl(columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData)
2234
2235
2236 end subroutine oop_Hto
2237
2238
2239 subroutine oop_HheightCoordObs()
2240 !
2241 !:Purpose: Compute simulated geometric-height based observations
2242 ! such as Radar Doppler velocity.
2243 ! It returns Hdx in OBS_WORK
2244 implicit none
2245
2246 ! Locals:
2247 integer :: bodyIndex, headerIndex, levelIndex, bufrCode, layerIndex, familyIndex
2248 integer :: bodyIndexStart, bodyIndexEnd, bodyIndex1
2249 real(8) :: levelAltLow, levelAltHigh, HDX, Azimuth, obsAltitude
2250 real(8) :: uuLow, uuHigh, vvLow, vvHigh, vInterpWeightHigh, vInterpWeightLow
2251 real(8) :: anlIncValueLow, anlIncValueHigh, trlValueLow, trlValueHigh
2252 real(8), pointer :: du_column(:), dv_column(:), height_column(:)
2253 integer, parameter :: NUMFAMILY=3
2254 character(len=2) :: listFamily(NUMFAMILY), cfam
2255
2256 listFamily(1) = 'RA' ! Doppler velocity (Radial Wind) burf_radvel
2257 listFamily(2) = 'PR' ! Dew point difference burf_nees
2258 listFamily(3) = 'AL' ! Aladin HLOS wind burf_neal
2259 ! Loop over all family
2260 FAMILY: do familyIndex=1, NUMFAMILY
2261
2262 ! Loop over all header indices
2263 call obs_set_current_header_list(obsSpaceData, listFamily(familyIndex))
2264 HEADER: do
2265 headerIndex = obs_getHeaderIndex(obsSpaceData)
2266 if ( headerIndex < 0 ) exit HEADER
2267
2268 if ( listFamily(familyIndex) == 'RA' ) then
2269
2270 ! Azimuth of the radar beam
2271 azimuth = obs_headElem_r(obsSpaceData, OBS_RZAM, headerIndex)
2272
2273 end if
2274
2275 ! Local vector state
2276 du_column => col_getColumn(columnAnlInc, headerIndex, 'UU')
2277 dv_column => col_getColumn(columnAnlInc, headerIndex, 'VV')
2278
2279 ! Loop over all body indices
2280 call obs_set_current_body_list(obsSpaceData, headerIndex)
2281 BODY: do
2282 bodyIndex = obs_getBodyIndex(obsSpaceData)
2283 if ( bodyIndex < 0 ) exit BODY
2284
2285 ! only process observations flagged to be assimilated
2286 if ( obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated ) cycle BODY
2287
2288 ! Check that this observation has the expected bufr element ID
2289 bufrCode = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex)
2290 if ( bufrCode == BUFR_logRadarPrecip ) cycle BODY
2291
2292 if ( bufrCode == bufr_radvel ) then
2293 ! Operator Hx for tangential wind is linear in x
2294 ! that is:
2295 ! h(x) = Hx
2296 ! and
2297 ! h(xb) + Hdx = Hxb + Hdx
2298 !
2299 ! H includes vertical interpolation
2300 ! and projection of U and V wind components along the direction of the beam
2301 !
2302 ! In matrix form for one Doppler velocity observation:
2303 ! VDoppler = Hx = [ iwHigh*sin(az) iwHigh*cos(az) iwLow*sin(az) iwLow*cos(az) ][ uuHigh ]
2304 ! [ vvHigh ]
2305 ! [ uuLow ]
2306 ! [ vvLow ]
2307 ! such that
2308 ! VDoppler = iwHigh*sin(az)*uuHigh
2309 ! + iwHigh*cos(az)*vvHigh
2310 ! + iwLow *sin(az)*uuLow
2311 ! + iwLow *cos(az)*vvLow
2312 !
2313 ! With
2314 ! az = beam azimuth (radians, met convention)
2315 ! iwHigh = Interpolation Weight for model level just above observation
2316 ! iwLow = " below "
2317 ! uuHigh and vvHigh = wind components on model level just above observation
2318 ! uuLow and vvLow = " below "
2319 !
2320 ! The dependence of model levels on surface pressure is neglected here
2321
2322 ! OBS_LYR returns the index of the model level just above the observation
2323 ! level=1 is the highest level such that level+1 is lower
2324 levelIndex = obs_bodyElem_i(obsSpaceData, OBS_LYR, bodyIndex)
2325 levelAltHigh = col_getHeight(columnTrlOnAnlIncLev, levelIndex, headerIndex, 'MM')
2326 levelAltLow = col_getHeight(columnTrlOnAnlIncLev, levelIndex+1, headerIndex, 'MM')
2327 obsAltitude = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
2328
2329 !vertical interpolation weights
2330 vInterpWeightHigh = (obsAltitude - levelAltLow)/(levelAltHigh - levelAltLow)
2331 vInterpWeightLow = 1.0D0 - vInterpWeightHigh
2332
2333 HDX = vInterpWeightHigh*sin(azimuth)*du_column(levelIndex) &
2334 + vInterpWeightHigh*cos(azimuth)*dv_column(levelIndex) &
2335 + vInterpWeightLow *sin(azimuth)*du_column(levelIndex+1) &
2336 + vInterpWeightLow *cos(azimuth)*dv_column(levelIndex+1)
2337
2338 ! Store HDX in OBS_WORK
2339 call obs_bodySet_r(obsSpaceData, OBS_WORK, bodyIndex, HDX)
2340
2341 else
2342
2343 cfam = obs_getfamily(obsSpaceData,headerIndex_opt=headerIndex)
2344 write(*,*) 'CANNOT ASSIMILATE OBSERVATION!!!', &
2345 'bufrCode =', bufrCode, 'cfam =', cfam
2346 call utl_abort('oop_HheightCoordObs')
2347
2348 end if
2349 end do BODY
2350 end do HEADER
2351 end do FAMILY
2352
2353 end subroutine oop_HheightCoordObs
2354
2355 !--------------------------------------------------------------------------
2356
2357 subroutine oop_Hro( initializeLinearization_opt )
2358 !
2359 !:Purpose: Compute the tangent operator for GPSRO observations.
2360 !
2361 implicit none
2362
2363 ! Arguments:
2364 logical, optional, intent(in) :: initializeLinearization_opt
2365
2366 ! Locals:
2367 real(8) :: ZMHXL
2368 real(8) :: DX(gps_ncvmx)
2369 integer :: IDATYP
2370 integer :: JL, JV, NGPSLEV, NWNDLEV
2371 integer :: headerIndex, bodyIndex, iProfile
2372 logical :: ASSIM
2373 integer :: NH, NH1
2374
2375 ! Initializations
2376 NGPSLEV = col_getNumLev(columnAnlInc,'TH')
2377 NWNDLEV = col_getNumLev(columnAnlInc,'MM')
2378
2379 ! call to calculate the GPSRO Jacobians
2380 call oop_calcGPSROJacobian( columnTrlOnAnlIncLev, obsSpaceData, &
2381 initializeLinearization_opt=initializeLinearization_opt )
2382
2383 ! Loop over all header indices of the 'RO' family (Radio Occultation)
2384 ! Set the header list (start at the beginning of the list)
2385 call obs_set_current_header_list(obsSpaceData,'RO')
2386 HEADER: do
2387 headerIndex = obs_getHeaderIndex(obsSpaceData)
2388 if (headerIndex < 0) exit HEADER
2389
2390 ! Process only refractivity data (codtyp 169)
2391 IDATYP = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
2392 DATYP: if ( IDATYP == 169 ) THEN
2393
2394 ! Scan for requested data values of the profile, and count them
2395 ASSIM = .FALSE.
2396 NH = 0
2397
2398 ! Loop over all body indices for this headerIndex:
2399 ! (start at the beginning of the list)
2400 call obs_set_current_body_list(obsSpaceData, headerIndex)
2401 BODY: do
2402 bodyIndex = obs_getBodyIndex(obsSpaceData)
2403 if (bodyIndex < 0) exit BODY
2404 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) THEN
2405 ASSIM = .TRUE.
2406 NH = NH + 1
2407 end if
2408 end do BODY
2409
2410 ! If assimilations are requested, prepare and apply the observation operator
2411 ASSIMILATE: if (ASSIM) THEN
2412 iProfile=gps_iprofile_from_index(headerIndex)
2413
2414 ! Local vector state
2415 do JL = 1, NGPSLEV
2416 DX ( JL) = col_getElem(columnAnlInc,JL,headerIndex,'TT')
2417 DX (NGPSLEV+JL) = col_getElem(columnAnlInc,JL,headerIndex,'HU')
2418 end do
2419 DX (2*NGPSLEV+1:3*NGPSLEV) = col_getColumn(columnAnlInc,headerIndex,'Z_T')
2420 DX (3*NGPSLEV+1:4*NGPSLEV) = col_getColumn(columnAnlInc,headerIndex,'P_T')
2421
2422 ! Perform the (H(xb)DX-Y') operation
2423 ! Loop over all body indices for this headerIndex:
2424 NH1 = 0
2425 call obs_set_current_body_list(obsSpaceData, headerIndex)
2426 BODY_3: do
2427 bodyIndex = obs_getBodyIndex(obsSpaceData)
2428 if (bodyIndex < 0) exit BODY_3
2429 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) THEN
2430 NH1 = NH1 + 1
2431
2432 ! Evaluate H(xb)DX
2433 ZMHXL = 0.d0
2434 do JV = 1, 4*NGPSLEV
2435 ZMHXL = ZMHXL + oop_vRO_Jacobian(iProfile,NH1,JV) * DX(JV)
2436 end do
2437
2438 ! Store in CMA
2439 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex, ZMHXL)
2440 end if
2441 end do BODY_3
2442 end if ASSIMILATE
2443 end if DATYP
2444 end do HEADER
2445
2446
2447 end subroutine oop_Hro
2448
2449 !--------------------------------------------------------------------------
2450
2451 subroutine oop_Hgp( initializeLinearization_opt )
2452 !
2453 !:Purpose: Compute H'dx for all GPS ZTD observations
2454 ! oop_Hgp TL of DOBSGPSGB (Jo for GB-GPS ZTD observations)
2455 implicit none
2456
2457 ! Arguments:
2458 logical, optional, intent(in) :: initializeLinearization_opt
2459
2460 ! Locals:
2461 real(8) :: ZHX
2462 real(8) :: DX(gps_ncvmx)
2463 integer :: headerIndex, bodyIndex
2464 integer :: JL, NFLEV, status, iztd, icount
2465 logical :: ASSIM
2466 character(len=12) :: cstnid
2467
2468 NFLEV = col_getNumLev(columnTrlOnAnlIncLev,'TH')
2469
2470 icount = 0
2471
2472 ! call to calculate the GPSGB Jacobians
2473 call oop_calcGPSGBJacobian( columnTrlOnAnlIncLev, obsSpaceData, &
2474 initializeLinearization_opt=initializeLinearization_opt )
2475
2476 ! loop over all header indices of the 'GP' family (GPS observations)
2477 call obs_set_current_header_list(obsSpaceData,'GP')
2478 HEADER: do
2479 headerIndex = obs_getHeaderIndex(obsSpaceData)
2480 if (headerIndex < 0) exit HEADER
2481
2482 ! Scan for ZTD assimilation at this location
2483 ASSIM = .FALSE.
2484 ! loop over all body indices for this headerIndex
2485 call obs_set_current_body_list(obsSpaceData, headerIndex)
2486 BODY: do
2487 bodyIndex = obs_getBodyIndex(obsSpaceData)
2488 if (bodyIndex < 0) exit BODY
2489
2490 if ( (obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex) == 189) &
2491 .and. (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == bufr_nezd) &
2492 .and. (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) ) then
2493 ASSIM = .TRUE.
2494 end if
2495 end do BODY
2496
2497 ! If ZTD assimilation, apply the TL observation operator
2498 if ( ASSIM ) THEN
2499 iztd = gps_iztd_from_index(headerIndex)
2500 if ( iztd < 1 .or. iztd > gps_gb_numZTD ) then
2501 call utl_abort('oop_Hgp: ERROR: index from gps_iztd_from_index() is out of range!')
2502 end if
2503
2504 ! Local vector state (analysis increments)
2505 do JL = 1, NFLEV
2506 DX (JL) = col_getElem(columnAnlInc,JL,headerIndex,'TT')
2507 DX (NFLEV+JL) = col_getElem(columnAnlInc,JL,headerIndex,'HU')
2508 end do
2509 DX (2*NFLEV+1:3*NFLEV) = col_getColumn(columnAnlInc,headerIndex,'Z_T')
2510 DX (3*NFLEV+1:4*NFLEV) = col_getColumn(columnAnlInc,headerIndex,'P_T')
2511
2512 ! Evaluate H'(xb)*dX
2513 ZHX = 0.D0
2514 do JL = 1, 4*NFLEV
2515 ZHX = ZHX + oop_vZTD_Jacobian(iztd,JL)*DX(JL)
2516 end do
2517
2518 ! Store ZHX = H'dx in OBS_WORK
2519 ! loop over all body indices for this headerIndex
2520 call obs_set_current_body_list(obsSpaceData, headerIndex)
2521 BODY_2: do
2522 bodyIndex = obs_getBodyIndex(obsSpaceData)
2523 if (bodyIndex < 0) exit BODY_2
2524
2525 if ( (obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex) == 189) &
2526 .and. (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == bufr_nezd) &
2527 .and. (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) ) then
2528 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex, ZHX)
2529 icount = icount + 1
2530 if ( icount <= 3 .and. gps_gb_LTESTOP ) then
2531 cstnid = obs_elem_c(obsSpaceData,'STID',headerIndex)
2532 write(*,*) iztd, cstnid
2533 write(*,*) 'JAC(ncv) = ', (oop_vZTD_Jacobian(iztd,JL),JL=1,4*NFLEV)
2534 write(*,*) 'DTT(JL) = ', (DX(JL),JL=1,NFLEV)
2535 write(*,*) 'DHU(JL) = ', (DX(JL),JL=NFLEV+1,2*NFLEV)
2536 write(*,*) 'DAL(JL) = ', (DX(JL),JL=2*NFLEV+1,3*NFLEV)
2537 write(*,*) 'DP (JL) = ', (DX(JL),JL=3*NFLEV+1,4*NFLEV)
2538 write(*,*) 'ZHX (mm) = ', ZHX*1000.D0
2539 end if
2540 end if
2541 end do BODY_2
2542
2543 end if ! ASSIM
2544
2545 end do HEADER
2546
2547 ! WRITE(*,*) 'oop_Hgp: Number of ZTD data locations with obs_bodySet_r(OBS_OMA) = ', icount
2548
2549 ! WRITE(*,*)'EXIT oop_Hgp'
2550
2551
2552 end subroutine oop_Hgp
2553
2554 !--------------------------------------------------------------------------
2555
2556 subroutine oop_Hchm()
2557 !:Purpose: Compute simulated chemical constituents observations from profiled model
2558 ! increments, and returns Hdx in OBS_WORK
2559 implicit none
2560
2561 if (.not.obs_famExist(obsSpaceData,'CH', localMPI_opt = .true. )) return
2562
2563 call oopc_CHobsoperators(columnTrlOnAnlIncLev,obsSpaceData,'tl',columnAnlInc_opt=columnAnlInc) ! 'tl' for tangent linear operator
2564
2565 end subroutine oop_Hchm
2566
2567 end subroutine oop_Htl
2568
2569 !--------------------------------------------------------------------------
2570 ! oop_Had
2571 !--------------------------------------------------------------------------
2572 subroutine oop_Had(columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData, initializeLinearization_opt)
2573 !
2574 !:Purpose: Call the several adjoint of observation operators
2575 implicit none
2576
2577 ! Arguments:
2578 type(struct_columnData), intent(inout) :: columnAnlInc
2579 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
2580 type(struct_obs) , intent(inout) :: obsSpaceData
2581 logical, optional , intent(in) :: initializeLinearization_opt
2582
2583 ! Locals:
2584 type(struct_vco), pointer :: vco_anl
2585 logical, save :: initializeLinearization = .true.
2586
2587 if ( mmpi_myid == 0 ) then
2588 write(*,*)'OOP_HT- Adjoint of linearized observation operators'
2589 end if
2590
2591 vco_anl => col_getVco(columnTrlOnAnlIncLev)
2592
2593 ! Re-linearlize if it is asked for
2594 if ( present(initializeLinearization_opt) ) then
2595 initializeLinearization = initializeLinearization_opt
2596 end if
2597
2598 ! Find interpolation layer in model profiles (used by several operators)
2599 if ( initializeLinearization ) then
2600 if ( col_getNumLev(columnTrlOnAnlIncLev,'MM') > 1 ) call oop_vobslyrs(columnTrlOnAnlIncLev, obsSpaceData, beSilent=.false.)
2601 initializeLinearization = .false.
2602 end if
2603
2604 call oop_HTchm
2605
2606 if ( gps_gb_numZTD > 0 ) then
2607 call oop_HTgp( initializeLinearization_opt=initializeLinearization_opt )
2608 end if
2609
2610 call oop_HTheighCoordObs()
2611
2612 call oop_HTro( initializeLinearization_opt=initializeLinearization_opt )
2613
2614 call oop_HTto
2615
2616 call oop_HTsf
2617
2618 call oop_HTpp
2619
2620 call oop_HTsst
2621
2622 call oop_HTice
2623
2624 call oop_HThydro()
2625
2626
2627 CONTAINS
2628
2629 !--------------------------------------------------------------------------
2630
2631 subroutine oop_HTpp
2632 !:Purpose: Adjoint of the "vertical" interpolation, based on vint3d,
2633 ! for "UPPER AIR" data files.
2634 implicit none
2635
2636 ! Locals:
2637 integer levIndexBot,levIndexTop,bufrCode
2638 real(8) ZRES
2639 real(8) ZWB,ZWT
2640 real(8) ZLEV,ZPT,ZPB,ZDADPS,ZPRESBPB,ZPRESBPT
2641 integer headerIndex,layerIndex,nlev_T
2642 integer bodyIndex,INDEX_FAMILY
2643 real(8) trlValueTop,trlValueBot
2644 real(8), pointer :: all_column(:),tt_column(:),hu_column(:),p_column(:)
2645 real(8) :: delPT,delPB
2646 integer, parameter :: numFamily=3
2647 character(len=2) :: list_family(numFamily)
2648 character(len=4) :: varLevel
2649
2650 list_family(1) = 'UA'
2651 list_family(2) = 'AI'
2652 list_family(3) = 'SW'
2653
2654 FAMILY: do index_family=1,numFamily
2655
2656 call obs_set_current_body_list(obsSpaceData,list_family(index_family))
2657 BODY: do
2658 bodyIndex = obs_getBodyIndex(obsSpaceData)
2659 if (bodyIndex < 0) exit BODY
2660
2661 ! Process all data within the domain of the model
2662 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated .and. &
2663 obs_bodyElem_i(obsSpaceData,OBS_XTR,bodyIndex) == 0 .and. &
2664 obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) == 2 ) then
2665
2666 headerIndex = obs_bodyElem_i(obsSpaceData,OBS_HIND,bodyIndex)
2667 ZRES = obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
2668 ZLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
2669 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
2670 varLevel = vnl_varLevelFromVarnum(bufrCode)
2671 layerIndex = obs_bodyElem_i(obsSpaceData,OBS_LYR,bodyIndex)
2672 levIndexTop = layerIndex + col_getOffsetFromVarno(columnTrlOnAnlIncLev,bufrCode)
2673 levIndexBot = levIndexTop+1
2674 ZPT = col_getPressure(columnTrlOnAnlIncLev,layerIndex,headerIndex,varLevel)
2675 ZPB = col_getPressure(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,varLevel)
2676 delPT= col_getPressure(columnAnlInc,layerIndex ,headerIndex,varLevel)
2677 delPB= col_getPressure(columnAnlInc,layerIndex+1,headerIndex,varLevel)
2678 ZWB = LOG(ZLEV/ZPT)/LOG(ZPB/ZPT)
2679 ZWT = 1.0D0 - ZWB
2680
2681 ZDADPS = ( LOG(ZLEV/ZPB)*delPT/ZPT - &
2682 LOG(ZLEV/ZPT)*delPB/ZPB )/ &
2683 LOG(ZPB/ZPT)**2
2684
2685 all_column => col_getColumn(columnAnlInc,headerIndex)
2686 tt_column => col_getColumn(columnAnlInc,headerIndex,'TT')
2687 hu_column => col_getColumn(columnAnlInc,headerIndex,'HU')
2688 if ( varLevel == 'TH' ) then
2689 p_column => col_getColumn(columnAnlInc,headerIndex,'P_T')
2690 else if ( varLevel == 'MM' ) then
2691 p_column => col_getColumn(columnAnlInc,headerIndex,'P_M')
2692 end if
2693
2694 if (bufrCode.eq.bufr_nees) then
2695 call hutoes_ad(hu_column(layerIndex+1), &
2696 tt_column(layerIndex+1), &
2697 p_column(layerIndex+1), &
2698 ZWB*ZRES, &
2699 col_getElem(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'HU'), &
2700 col_getPressure(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'TH'))
2701 call hutoes_ad(hu_column(layerIndex ), &
2702 tt_column(layerIndex ), &
2703 p_column(layerIndex), &
2704 ZWT*ZRES, &
2705 col_getElem(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'HU'), &
2706 col_getPressure(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'TH'))
2707 trlValueBot=hutoes(col_getElem(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'HU'), &
2708 col_getElem(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'TT'), &
2709 col_getPressure(columnTrlOnAnlIncLev,layerIndex+1,headerIndex,'TH'))
2710 trlValueTop=hutoes(col_getElem(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'HU'), &
2711 col_getElem(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'TT'), &
2712 col_getPressure(columnTrlOnAnlIncLev,layerIndex ,headerIndex,'TH'))
2713 else
2714 all_column(levIndexBot) = all_column(levIndexBot) + ZWB*ZRES
2715 all_column(levIndexTop) = all_column(levIndexTop) + ZWT*ZRES
2716 trlValueBot=col_getElem(columnTrlOnAnlIncLev,levIndexBot,headerIndex)
2717 trlValueTop=col_getElem(columnTrlOnAnlIncLev,levIndexTop,headerIndex)
2718 end if
2719 p_column(layerIndex ) = p_column(layerIndex ) + (trlValueBot - trlValueTop) * &
2720 (LOG(ZLEV/ZPB)/ZPT/LOG(ZPB/ZPT)**2) * ZRES
2721 p_column(layerIndex+1) = p_column(layerIndex+1) - (trlValueBot - trlValueTop) * &
2722 (LOG(ZLEV/ZPT)/ZPB/LOG(ZPB/ZPT)**2) * ZRES
2723 end if
2724
2725 end do BODY
2726
2727 end do FAMILY
2728
2729 end subroutine oop_HTpp
2730
2731 !--------------------------------------------------------------------------
2732
2733 subroutine oop_HTsf
2734 !:Purpose: based on surfc1dz to build the adjoint of the
2735 ! vertical interpolation for SURFACE data files.
2736 implicit none
2737
2738 ! Locals:
2739 integer :: levIndexBot,levIndexTop
2740 real(8) :: ZRES
2741 real(8) :: ZWB, ZWT,coeffA,coeffB,ZATV,trlVirtTemp
2742 real(8) :: obsHeight, ZPT, ZPB, ZDADPS, ZDELPS, ZDELTV, deltaT
2743 real(8) :: trlValueBot
2744 real(8) :: trlUwind, trlVwind, sumSquare, trlWindSpeed, anlIncWindSpeed
2745 integer :: headerIndex, layerIndex, nlev, nlev_T
2746 integer :: bodyIndex, bufrCode, INDEX_FAMILY
2747 real(8), pointer :: all_column(:), tt_column(:), hu_column(:), ps_column(:), p_column(:)
2748 real(8), pointer :: du_column(:), dv_column(:)
2749 real(8) :: dPdPsfc
2750 integer, parameter :: numFamily=5
2751 character(len=2) :: list_family(numFamily)
2752 character(len=4) :: varLevel
2753
2754 !- Temperature lapse rate for extrapolation of height below model surface
2755 coeffB = 1.0d0/(MPC_RGAS_DRY_AIR_R8*temperatureLapseRate/ec_rg)
2756
2757 !
2758 !- 1. Fill in COMMVO by using the adjoint of the "vertical" interpolation
2759 !
2760 list_family(1) = 'UA'
2761 list_family(2) = 'SF'
2762 list_family(3) = 'SC'
2763 list_family(4) = 'GP'
2764 list_family(5) = 'RA'
2765
2766 FAMILY: do index_family=1,numFamily
2767
2768 call obs_set_current_body_list(obsSpaceData, list_family(index_family))
2769 BODY: do
2770 bodyIndex = obs_getBodyIndex(obsSpaceData)
2771 if (bodyIndex < 0) exit BODY
2772
2773 ! Process all data within the domain of the model
2774 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
2775 if ( bufrCode == bufr_nezd .or. bufrCode == bufr_radvel ) cycle BODY
2776 if ( (obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex) == 1) &
2777 .and. (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) &
2778 .and. (bufrCode == bufr_nets .or. bufrCode == bufr_neps &
2779 .or. bufrCode == bufr_nepn .or. bufrCode == bufr_ness &
2780 .or. bufrCode == bufr_neus .or. bufrCode == bufr_nevs &
2781 .or. bufrCode == bufr_vis .or. bufrCode == bufr_logVis &
2782 .or. bufrCode == bufr_gust .or. bufrCode == bufr_nefs &
2783 .or. bufrCode == bufr_radarPrecip .or. bufrCode == bufr_logRadarPrecip &
2784 .or. obs_bodyElem_i(obsSpaceData,OBS_XTR,bodyIndex) == 0) ) then
2785
2786 varLevel = vnl_varLevelFromVarnum(bufrCode)
2787 nlev = col_getNumLev(columnAnlInc,varLevel)
2788 nlev_T = col_getNumLev(columnAnlInc,'TH')
2789 headerIndex = obs_bodyElem_i(obsSpaceData,OBS_HIND,bodyIndex)
2790 bufrCode = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
2791 layerIndex = obs_bodyElem_i(obsSpaceData,OBS_LYR,bodyIndex)
2792 obsHeight = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
2793 levIndexTop = nlev - 1 + col_getOffsetFromVarno(columnTrlOnAnlIncLev,bufrCode)
2794 levIndexBot = levIndexTop+1
2795 ZRES = obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
2796
2797 if (bufrCode == bufr_nets .or. bufrCode == bufr_ness .or. &
2798 bufrCode == bufr_neus .or. bufrCode == bufr_nevs .or. &
2799 bufrCode == bufr_vis .or. bufrCode == bufr_logVis .or. &
2800 bufrCode == bufr_gust .or. &
2801 bufrCode == bufr_radarPrecip .or. bufrCode == bufr_logRadarPrecip) then
2802 if ( bufrCode == bufr_ness ) then
2803 tt_column => col_getColumn(columnAnlInc,headerIndex,'TT')
2804 hu_column => col_getColumn(columnAnlInc,headerIndex,'HU')
2805 p_column => col_getColumn(columnAnlInc,headerIndex,'P_T')
2806 call hutoes_ad(hu_column(nlev_T), &
2807 tt_column(nlev_T), &
2808 p_column(nlev_T), &
2809 ZRES, &
2810 col_getElem(columnTrlOnAnlIncLev,nlev_T,headerIndex,'HU'), &
2811 col_getPressure(columnTrlOnAnlIncLev,nlev_T,headerIndex,'TH'))
2812 else
2813 all_column => col_getColumn(columnAnlInc,headerIndex)
2814 all_column(levIndexBot) = all_column(levIndexBot) + ZRES
2815 end if
2816 else if ( bufrCode == bufr_nefs ) then
2817 du_column => col_getColumn(columnAnlInc,headerIndex,'UU')
2818 dv_column => col_getColumn(columnAnlInc,headerIndex,'VV')
2819 anlIncWindSpeed=obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
2820 sumSquare = col_getElem(columnTrlOnAnlIncLev,nlev,headerIndex,'UU')**2 + &
2821 col_getElem(columnTrlOnAnlIncLev,nlev,headerIndex,'VV')**2
2822 if ( sumSquare > 1.0d-10 ) then
2823 trlWindSpeed=sqrt(sumSquare)
2824 du_column(nlev) = du_column(nlev) + &
2825 anlIncWindSpeed*col_getElem(columnTrlOnAnlIncLev,nlev,headerIndex,'UU')/trlWindSpeed
2826 dv_column(nlev) = dv_column(nlev) + &
2827 anlIncWindSpeed*col_getElem(columnTrlOnAnlIncLev,nlev,headerIndex,'VV')/trlWindSpeed
2828 else
2829 trlWindSpeed=0.
2830 end if
2831 else if ( bufrCode == bufr_neps .or. bufrCode == bufr_nepn ) then
2832 tt_column => col_getColumn(columnAnlInc,headerIndex,'TT')
2833 hu_column => col_getColumn(columnAnlInc,headerIndex,'HU')
2834 ps_column => col_getColumn(columnAnlInc,headerIndex,'P0')
2835 trlVirtTemp = columnTrlOnAnlIncLev%OLTV(1,nlev_T,headerIndex)*col_getElem(columnTrlOnAnlIncLev,nlev_T,headerIndex,'TT')
2836 deltaT = temperatureLapseRate*(obsHeight-col_getHeight(columnTrlOnAnlIncLev,0,headerIndex,'SF'))
2837 coeffA = ((trlVirtTemp-deltaT)/trlVirtTemp)
2838 ZDELTV = (col_getElem(columnTrlOnAnlIncLev,1,headerIndex,'P0')*coeffB*coeffA**(coeffB-1)) &
2839 *(deltaT/(trlVirtTemp*trlVirtTemp))
2840 ZDELPS = coeffA**coeffB
2841 ZATV = ZDELTV*ZRES
2842 ps_column(1)= ps_column(1) + ZDELPS*ZRES
2843 tt_column(nlev_T) = tt_column(nlev_T) + &
2844 columnTrlOnAnlIncLev%OLTV(1,nlev_T,headerIndex)*ZATV
2845 hu_column(nlev_T) = hu_column(nlev_T) + &
2846 columnTrlOnAnlIncLev%OLTV(2,nlev_T,headerIndex)*ZATV
2847 else
2848
2849 call utl_abort('oop_HTsf: You have entered the twilight zone!')
2850
2851 end if
2852
2853 end if
2854
2855 end do BODY
2856
2857 end do FAMILY
2858
2859 end subroutine oop_HTsf
2860
2861 !--------------------------------------------------------------------------
2862
2863 subroutine oop_HTsst
2864 !:Purpose: Adjoint of the "vertical" interpolation for SST data
2865 implicit none
2866
2867 ! Locals:
2868 real(8) :: residual
2869 integer :: headerIndex, bodyIndex, bufrCode
2870 real(8), pointer :: columnTG(:)
2871 character(len=4) :: varName
2872
2873 call obs_set_current_body_list( obsSpaceData, 'TM' )
2874
2875 !$OMP PARALLEL DO PRIVATE(bodyIndex, bufrCode, varName, headerIndex, residual, columnTG)
2876 BODY: do bodyIndex = 1, obs_numBody( obsSpaceData )
2877 if ( obs_getFamily(obsSpaceData, bodyIndex_opt=bodyIndex) /= 'TM' ) cycle BODY
2878
2879 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
2880 if ( bufrCode /= bufr_sst ) cycle BODY
2881
2882 if ( col_varExist(columnAnlInc,'TM') ) then
2883 varName = 'TM'
2884 else
2885 varName = 'TG'
2886 end if
2887
2888 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated ) then
2889
2890 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
2891 residual = obs_bodyElem_r( obsSpaceData, OBS_WORK, bodyIndex )
2892 columnTG => col_getColumn(columnAnlInc, headerIndex, varName_opt = varName )
2893 columnTG(1) = columnTG(1) + residual
2894
2895 end if
2896
2897 end do BODY
2898 !$OMP END PARALLEL DO
2899
2900 end subroutine oop_HTsst
2901
2902 !--------------------------------------------------------------------------
2903
2904 subroutine oop_HThydro
2905 !:Purpose: Adjoint of the "vertical" interpolation for Hydrological data
2906 implicit none
2907
2908 ! Locals:
2909 real(8) :: residual
2910 integer :: headerIndex, bodyIndex, bufrCode
2911 real(8), pointer :: columnHY(:)
2912 character(len=4) :: varName
2913
2914 call obs_set_current_body_list( obsSpaceData, 'HY' )
2915
2916 BODY: do
2917
2918 bodyIndex = obs_getBodyIndex( obsSpaceData )
2919 if (bodyIndex < 0) exit BODY
2920
2921 ! Process all data within the domain of the model
2922 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
2923
2924 if ( bufrCode /= bufr_riverFlow ) cycle BODY
2925
2926 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated ) then
2927
2928 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
2929 residual = obs_bodyElem_r( obsSpaceData, OBS_WORK, bodyIndex )
2930 varName = vnl_varNameFromVarNum(bufrCode)
2931 columnHY => col_getColumn(columnAnlInc, headerIndex, varName_opt = varName )
2932 columnHY(1) = columnHY(1) + residual
2933 end if
2934
2935 end do BODY
2936
2937 end subroutine oop_HThydro
2938
2939 !--------------------------------------------------------------------------
2940
2941 subroutine oop_HTice
2942 !:Purpose: Adjoint of the "vertical" interpolation for ICE data
2943 implicit none
2944
2945 ! Locals:
2946 real(8) :: residual, scaling
2947 integer :: headerIndex, bodyIndex, bufrCode
2948 real(8), pointer :: columnGL(:)
2949 character(len=4) :: varName
2950
2951 call obs_set_current_body_list( obsSpaceData, 'GL' )
2952
2953 BODY: do
2954
2955 bodyIndex = obs_getBodyIndex( obsSpaceData )
2956 if (bodyIndex < 0) exit BODY
2957
2958 ! Process all data within the domain of the model
2959
2960 scaling = oop_iceScaling(obsSpaceData, bodyIndex)
2961
2962 if (scaling == 0.0d0) cycle BODY
2963
2964 if ( obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated ) then
2965 residual = scaling*obs_bodyElem_r( obsSpaceData, OBS_WORK, bodyIndex )
2966 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
2967 varName = vnl_varNameFromVarNum(bufrCode)
2968 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
2969 columnGL => col_getColumn( columnAnlInc, headerIndex, varName_opt = varName )
2970 columnGL(1) = columnGL(1) + residual
2971 end if
2972
2973 end do BODY
2974
2975 end subroutine oop_HTice
2976
2977 !--------------------------------------------------------------------------
2978
2979 subroutine oop_HTto
2980 !:Purpose: Adjoint of computation of residuals to the tovs observations
2981 implicit none
2982
2983 ! Locals:
2984 integer :: datestamp
2985
2986 if (.not.obs_famExist(obsSpaceData,'TO', localMPI_opt = .true. )) return
2987
2988 ! Adjoint of computing radiance
2989 datestamp = tim_getDatestamp()
2990 call tvs_fillProfiles(columnTrlOnAnlIncLev,obsSpaceData,datestamp,"tlad",.false.)
2991
2992 call tvslin_rttov_ad(columnAnlInc,columnTrlOnAnlIncLev,obsSpaceData)
2993
2994 end subroutine oop_HTto
2995
2996 !--------------------------------------------------------------------------
2997
2998 subroutine oop_HTro ( initializeLinearization_opt )
2999 !
3000 !:Purpose: Compute the adjoint operator for GPSRO observations.
3001 implicit none
3002
3003 ! Arguments:
3004 logical, optional, intent(in) :: initializeLinearization_opt
3005
3006 ! Locals:
3007 real(8) :: DPJO0(gps_ncvmx)
3008 real(8) :: DPJO1(gps_ncvmx)
3009 real(8) :: ZINC
3010 real(8), pointer :: tt_column(:),hu_column(:),height_column(:),p_column(:)
3011 integer :: IDATYP
3012 integer :: JL, NGPSLEV
3013 integer :: headerIndex, bodyIndex, iProfile
3014 logical :: ASSIM, LUSE
3015 integer :: NH, NH1
3016
3017 ! Initializations
3018 NGPSLEV=col_getNumLev(columnAnlInc,'TH')
3019
3020 ! call to calculate the GPSRO Jacobians
3021 call oop_calcGPSROJacobian( columnTrlOnAnlIncLev, obsSpaceData, &
3022 initializeLinearization_opt=initializeLinearization_opt )
3023
3024 ! Loop over all header indices of the 'RO' family (Radio Occultation)
3025 ! Set the header list (start at the beginning of the list)
3026 call obs_set_current_header_list(obsSpaceData,'RO')
3027 HEADER: do
3028 headerIndex = obs_getHeaderIndex(obsSpaceData)
3029 if (headerIndex < 0) exit HEADER
3030
3031 DPJO0 = 0.d0
3032
3033 ! Process only refractivity data (codtyp 169)
3034 IDATYP = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
3035 DATYP: if ( IDATYP == 169 ) THEN
3036
3037 ! Scan for requested data values of the profile, and count them
3038 ASSIM = .FALSE.
3039 NH = 0
3040
3041 ! Loop over all body indices for this headerIndex:
3042 ! (start at the beginning of the list)
3043 call obs_set_current_body_list(obsSpaceData, headerIndex)
3044 BODY: do
3045 bodyIndex = obs_getBodyIndex(obsSpaceData)
3046 if (bodyIndex < 0) exit BODY
3047
3048 LUSE=( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated )
3049 if ( LUSE ) THEN
3050 ASSIM = .TRUE.
3051 NH = NH + 1
3052 end if
3053 end do BODY
3054
3055 ! If assimilations are requested, prepare and apply the observation operator
3056 ASSIMILATE: if (ASSIM) THEN
3057 iProfile=gps_iprofile_from_index(headerIndex)
3058
3059 ! Perform the (H(xb)DX-Y')/S operation
3060 NH1 = 0
3061
3062 ! Loop over all body indices for this headerIndex:
3063 ! (start at the beginning of the list)
3064 call obs_set_current_body_list(obsSpaceData, headerIndex)
3065 BODY_3: do
3066 bodyIndex = obs_getBodyIndex(obsSpaceData)
3067 if (bodyIndex < 0) exit BODY_3
3068
3069 LUSE=( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated )
3070 if ( LUSE ) THEN
3071 NH1 = NH1 + 1
3072
3073 ! Normalized increment
3074 ZINC = obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
3075
3076 ! O-F Tested criteria:
3077 DPJO1(1:4*NGPSLEV) = ZINC * oop_vRO_Jacobian(iProfile,NH1,1:4*NGPSLEV)
3078
3079 ! Accumulate the gradient of the observation cost function:
3080 DPJO0(1:4*NGPSLEV) = DPJO0(1:4*NGPSLEV) + DPJO1(1:4*NGPSLEV)
3081 end if
3082 end do BODY_3
3083 end if ASSIMILATE
3084 end if DATYP
3085
3086 ! Store H* (HX - Z)/SIGMA in COMMVO
3087 tt_column => col_getColumn(columnAnlInc,headerIndex,'TT')
3088 hu_column => col_getColumn(columnAnlInc,headerIndex,'HU')
3089 height_column => col_getColumn(columnAnlInc,headerIndex,'Z_T')
3090 p_column => col_getColumn(columnAnlInc,headerIndex,'P_T')
3091 do JL = 1, NGPSLEV
3092 tt_column(JL) = DPJO0(JL)
3093 hu_column(JL) = DPJO0(JL+NGPSLEV)
3094 height_column(JL) = DPJO0(JL+2*NGPSLEV)
3095 p_column(JL) = DPJO0(JL+3*NGPSLEV)
3096 end do
3097 end do HEADER
3098
3099 end subroutine oop_HTro
3100
3101
3102 !--------------------------------------------------------------------------
3103
3104 subroutine oop_HTgp( initializeLinearization_opt )
3105 !
3106 !:Purpose: Compute Ht*grad(Jo) for all GPS ZTD observations
3107 !
3108 !:Note: ZTD Jacobians are computed and stored in oop_Hgp (first iter.)
3109 implicit none
3110
3111 ! Arguments:
3112 logical, optional, intent(in) :: initializeLinearization_opt
3113
3114 ! Locals:
3115 real(8) :: DPJO0(gps_ncvmx)
3116 real(8) :: JAC(gps_ncvmx)
3117 real(8) :: ZINC
3118 integer :: JL, NFLEV, iztd
3119 integer :: headerIndex, bodyIndex, icount
3120 logical :: ASSIM
3121 real(8), pointer :: tt_column(:),hu_column(:),height_column(:),p_column(:)
3122
3123 NFLEV = col_getNumLev(columnTrlOnAnlIncLev,'TH')
3124
3125 ! call to calculate the GPSGB Jacobians
3126 call oop_calcGPSGBJacobian( columnTrlOnAnlIncLev, obsSpaceData, &
3127 initializeLinearization_opt=initializeLinearization_opt )
3128
3129 ! loop over all header indices of the 'GP' family (GPS observations)
3130 ! Set the header list & start at the beginning of the list
3131 call obs_set_current_header_list(obsSpaceData,'GP')
3132
3133 icount = 0
3134
3135 HEADER: do
3136 headerIndex = obs_getHeaderIndex(obsSpaceData)
3137 if (headerIndex < 0) exit HEADER
3138
3139 DPJO0(:) = 0.0D0
3140 JAC(:) = 0.0D0
3141
3142 ! Scan for requested ZTD assimilation
3143 ASSIM = .FALSE.
3144
3145 ! loop over all body indices (still in the 'GP' family)
3146 ! Set the body list & start at the beginning of the list)
3147 call obs_set_current_body_list(obsSpaceData, headerIndex)
3148 BODY: do
3149 bodyIndex = obs_getBodyIndex(obsSpaceData)
3150 if (bodyIndex < 0) exit BODY
3151
3152 if ( (obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex) == 189) &
3153 .and. (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == bufr_nezd) &
3154 .and. (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) ) then
3155 ASSIM = .TRUE.
3156 end if
3157 end do BODY
3158 !
3159 if (ASSIM) THEN
3160
3161 icount = icount + 1
3162 iztd = gps_iztd_from_index(headerIndex)
3163 if ( iztd < 1 .or. iztd > gps_gb_numZTD ) then
3164 call utl_abort('oop_HTgp: ERROR: index from gps_iztd_from_index() is out of range!')
3165 end if
3166
3167 do JL = 1, 4*NFLEV
3168 JAC(JL) = oop_vZTD_Jacobian(iztd,JL)
3169 end do
3170
3171 ! Get Ht*grad(HeaderIndex) = Ht*(H'dx - d)/sigma_o^2
3172 ! loop over all body indices (still in the 'GP' family)
3173 ! Start at the beginning of the list)
3174 call obs_set_current_body_list(obsSpaceData, headerIndex)
3175 BODY_2: do
3176 bodyIndex = obs_getBodyIndex(obsSpaceData)
3177 if (bodyIndex < 0) exit BODY_2
3178 if ( (obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex ) == 189) &
3179 .and. (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == bufr_nezd) &
3180 .and. (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) ) then
3181 ZINC = obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
3182
3183 ! Accumulate the gradient of the observation cost function
3184 DPJO0(1:4*NFLEV) = ZINC * oop_vZTD_Jacobian(iztd,:)
3185 end if
3186 end do BODY_2
3187
3188 ! Store Ht*grad(HeaderIndex) in COMMVO
3189 tt_column => col_getColumn(columnAnlInc,headerIndex,'TT')
3190 hu_column => col_getColumn(columnAnlInc,headerIndex,'HU')
3191 height_column => col_getColumn(columnAnlInc,headerIndex,'Z_T')
3192 p_column => col_getColumn(columnAnlInc,headerIndex,'P_T')
3193 do JL = 1, NFLEV
3194 tt_column(JL) = DPJO0(JL)
3195 hu_column(JL) = DPJO0(JL+NFLEV)
3196 height_column(JL) = DPJO0(JL+2*NFLEV)
3197 p_column (JL) = DPJO0(JL+3*NFLEV)
3198 end do
3199
3200 end if ! ASSIM
3201
3202 end do HEADER
3203
3204 end subroutine oop_HTgp
3205
3206 ! --------------------------------------------------------------------------
3207
3208 subroutine oop_HTheighCoordObs()
3209 !:Purpose: Compute the adjoint operator of the
3210 ! vertical interpolation of geometric-height based data.
3211 ! Including Radar Doppler velocity data.
3212 !
3213 ! delta x is updated by this routine
3214 implicit none
3215
3216 ! Locals:
3217 integer :: bodyIndex, headerIndex, levelIndex, bufrCode, familyIndex
3218 integer :: bodyIndexStart, bodyIndexEnd, bodyIndex1
3219 real(8) :: levelAltLow, levelAltHigh, azimuth, obsAltitude, HDX
3220 real(8) :: anlIncValueLow, anlIncValueHigh, trlValueLow, trlValueHigh
3221 real(8) :: vInterpWeightLow, vInterpWeightHigh
3222 real(8), pointer :: du_column(:), dv_column(:), height_column(:)
3223 integer, parameter :: NUMFAMILY=3
3224 character(len=2) :: listFamily(NUMFAMILY), cfam
3225
3226 listFamily(1) = 'RA' ! Doppler velocity (Radial Wind) burf_radvel
3227 listFamily(2) = 'PR' ! Dew point difference burf_nees
3228 listFamily(3) = 'AL' ! Aladin HLOS wind burf_neal
3229
3230 ! Loop over all families
3231 FAMILY: do familyIndex = 1, NUMFAMILY
3232
3233 ! Loop over header indices
3234 call obs_set_current_header_list(obsSpaceData, listFamily(familyIndex))
3235 HEADER: do
3236 headerIndex = obs_getHeaderIndex(obsSpaceData)
3237 if ( headerIndex < 0 ) exit HEADER
3238
3239 if ( listFamily(familyIndex) == 'RA' ) then
3240 ! Azimuth of the radar beam
3241 azimuth = obs_headElem_r(obsSpaceData, OBS_RZAM, headerIndex)
3242 end if
3243
3244 ! Local vector state
3245 du_column => col_getColumn(columnAnlInc, headerIndex, 'UU')
3246 dv_column => col_getColumn(columnAnlInc, headerIndex, 'VV')
3247
3248 ! Loop over body indices
3249 call obs_set_current_body_list(obsSpaceData, headerIndex)
3250 BODY: do
3251 bodyIndex = obs_getBodyIndex(obsSpaceData)
3252 if ( bodyIndex < 0 ) exit BODY
3253 ! only process observations flagged to be assimilated
3254 if ( obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) /= obs_assimilated ) cycle BODY
3255
3256 ! Check that this observation has the expected bufr element ID
3257 bufrCode = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex)
3258 if ( bufrCode == BUFR_logRadarPrecip ) cycle BODY
3259
3260 if ( bufrCode == bufr_radvel ) then
3261 ! OBS_LYR returns the index of the model level just above the observation
3262 ! level=1 is the highest level such that level+1 is lower
3263 levelIndex = obs_bodyElem_i(obsSpaceData, OBS_LYR, bodyIndex)
3264 levelAltHigh = col_getHeight(columnTrlOnAnlIncLev, levelIndex, headerIndex, 'MM')
3265 levelAltLow = col_getHeight(columnTrlOnAnlIncLev, levelIndex+1, headerIndex, 'MM')
3266 obsAltitude = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
3267
3268 ! HDX from store
3269 HDX = obs_bodyElem_r(obsSpaceData, OBS_WORK, bodyIndex)
3270 vInterpWeightHigh = (obsAltitude - levelAltLow)/(levelAltHigh - levelAltLow)
3271 vInterpWeightLow = 1.0D0 - vInterpWeightHigh
3272
3273 ! see oop_HheightCoordObs for the description of the H operator for radar Doppler velocity
3274 !
3275 ! In matrix form, adjoint of H operator for one observation:
3276 !
3277 ! delta x = Ht TangentialWindResidual
3278 !
3279 ! [ deltaUUHigh ] = [ iwHigh*sin(az) ][ residual ]
3280 ! [ deltaVVigh ] [ iwHigh*cos(az) ]
3281 ! [ deltaUULow ] [ iwLow *sin(az) ]
3282 ! [ deltaVVLow ] [ iwLow *cos(az) ]
3283 !
3284 ! With
3285 ! az = beam azimuth (radians, met convention)
3286 ! iwHigh = Interpolation Weight for model level just above observation
3287 ! iwLow = " below "
3288 ! residual = OmP - Hdx
3289 ! deltaUUHigh and deltaVVHigh = wind components on model level just above observation
3290 ! deltaUULow and deltaVVLow = " below "
3291
3292 !delta x is updated
3293 du_column(levelIndex ) = du_column(levelIndex ) + vInterpWeightHigh*HDX*sin(azimuth)
3294 dv_column(levelIndex ) = dv_column(levelIndex ) + vInterpWeightHigh*HDX*cos(azimuth)
3295 du_column(levelIndex+1) = du_column(levelIndex+1) + vInterpWeightLow *HDX*sin(azimuth)
3296 dv_column(levelIndex+1) = dv_column(levelIndex+1) + vInterpWeightLow *HDX*cos(azimuth)
3297
3298 else
3299
3300 cfam = obs_getfamily(obsSpaceData,headerIndex_opt=headerIndex)
3301 write(*,*) 'CANNOT ASSIMILATE OBSERVATION!!!', &
3302 'bufrCode =', bufrCode, 'cfam =', cfam
3303 call utl_abort('oop_HTheighCoordObs')
3304
3305 end if
3306 end do BODY
3307 end do HEADER
3308 end do FAMILY
3309
3310 end subroutine oop_HTheighCoordObs
3311
3312 !--------------------------------------------------------------------------
3313
3314 subroutine oop_HTchm
3315 !:Purpose: Compute H^T * R^-1 (OmP-Hdx) for all CH observations
3316 implicit none
3317
3318 if (.not.obs_famExist(obsSpaceData,'CH', localMPI_opt = .true. )) return
3319
3320 call oopc_CHobsoperators(columnTrlOnAnlIncLev,obsSpaceData,'adjoint', &
3321 columnAnlInc_opt=columnAnlInc) ! 'adjoint' for adjoint of the tangent linear operator
3322
3323 end subroutine oop_HTchm
3324
3325 end subroutine oop_Had
3326
3327 !--------------------------------------------------------------------------
3328 ! HUtoES
3329 !--------------------------------------------------------------------------
3330 function HUtoES( hu, tt, pressure ) result( es )
3331 !:Purpose:
3332 ! to calculate the dew point depression from specific
3333 ! humidity, temperature and pressure. No ice phase
3334 ! is permitted and the pressure vector is given.
3335 implicit none
3336
3337 ! Arguments:
3338 real(8), intent(in) :: hu
3339 real(8), intent(in) :: tt
3340 real(8), intent(in) :: pressure
3341 ! Result:
3342 real(8) :: es
3343
3344 ! Locals:
3345 real(8) :: husat, td
3346
3347 ! get the saturated vapor pressure from specific humidity
3348 husat = phf_foefq8(hu,pressure)
3349
3350 ! now the dewpoint temperature
3351 td = phf_fotw8(husat)
3352
3353 ! finally the dewpoint depression
3354 es = min(tt-td,MPC_MAXIMUM_ES_R8)
3355
3356 end function HUtoES
3357
3358 !--------------------------------------------------------------------------
3359 ! HUtoES_tl
3360 !--------------------------------------------------------------------------
3361 function HUtoES_tl(HU_inc,TT_inc,P_inc,HU_trl,PRES_trl) result(ES_inc)
3362 !:Purpose: TLM VERSION
3363 ! to calculate the dew point depression from specific
3364 ! humidity, temperature and pressure. No ice phase
3365 ! is permitted and the pressure vector is given.
3366 implicit none
3367
3368 ! Arguments:
3369 real(8), intent(in) :: HU_inc
3370 real(8), intent(in) :: TT_inc
3371 real(8), intent(in) :: P_inc
3372 real(8), intent(in) :: HU_trl
3373 real(8), intent(in) :: PRES_trl
3374 ! Result:
3375 real(8) :: ES_inc
3376
3377 ! Locals:
3378 real(8) :: ZE, ZTD, dTDdE, ZQBRANCH
3379 real(8) :: dESdLQ, dESdTT, dESdP
3380
3381 dESdTT = 1.0d0
3382
3383 !- Forward calculations of saturation vapour pressure and dewpoint temperature
3384 ! and adjoint of vapour pressure from adjoint of dewpoint temperature
3385 ZE = phf_FOEFQ8(HU_trl, PRES_trl)
3386 ZTD = phf_FOTW8 (ZE)
3387 dTDdE= phf_FODTW8(ZTD,ZE)
3388
3389 !- adjoint of temp. specific humidity and surface pressure due to changes in vapour pressure
3390 ZQBRANCH = phf_FQBRANCH(HU_trl)
3391
3392 dESdLQ = - ZQBRANCH*phf_FOEFQA(1.0d0,dTDdE,HU_trl,PRES_trl)
3393
3394 dESdP = - ZQBRANCH*phf_FOEFQPSA(1.0d0,dTDdE,HU_trl,1.0d0)- &
3395 (1.D0-ZQBRANCH)*(dTDdE*1.0d0)
3396
3397 ES_inc = dESdLQ*HU_inc/HU_trl + dESdP*P_inc + dESdTT*TT_inc
3398
3399 end function HUtoES_tl
3400
3401 !--------------------------------------------------------------------------
3402 ! HUtoES_ad
3403 !--------------------------------------------------------------------------
3404 subroutine HUtoES_ad(HU_inc,TT_inc,P_inc,ES_inc,HU_trl,PRES_trl)
3405 ! Purpose: ADJOINT VERSION
3406 ! to calculate the dew point depression from specific
3407 ! humidity, temperature and pressure. No ice phase
3408 ! is permitted and the pressure vector is given.
3409 implicit none
3410
3411 ! Arguments:
3412 real(8), intent(inout) :: HU_inc
3413 real(8), intent(inout) :: TT_inc
3414 real(8), intent(inout) :: P_inc
3415 real(8), intent(in) :: ES_inc
3416 real(8), intent(in) :: HU_trl
3417 real(8), intent(in) :: PRES_trl
3418
3419 ! Locals:
3420 real(8) :: ZE,ZTD,dTDdE,ZQBRANCH
3421 real(8) :: dESdLQ,dESdTT,dESdP
3422
3423 dESdTT = 1.0d0
3424
3425 !- Forward calculations of saturation vapour pressure and dewpoint temperature
3426 ! and adjoint of vapour pressure from adjoint of dewpoint temperature
3427 ZE = phf_FOEFQ8(HU_trl, PRES_trl)
3428
3429 ZTD=phf_FOTW8(ZE)
3430 dTDdE=phf_FODTW8(ZTD,ZE)
3431
3432 !- adjoint of temp. specific humidity and surface pressure due to changes in vapour pressure
3433 ZQBRANCH = phf_FQBRANCH(HU_trl)
3434 dESdLQ = - ZQBRANCH*phf_FOEFQA(1.0d0,dTDdE,HU_trl,PRES_trl)
3435
3436 dESdP = - ZQBRANCH*phf_FOEFQPSA(1.0d0,dTDdE,HU_trl,1.0d0)- &
3437 (1.D0-ZQBRANCH)*(dTDdE*1.0d0)
3438
3439 ! TLM: ES_inc = dESdLQ*HU_inc/HU_trl + dESdP*P_inc + dESdTT*TT_inc
3440 ! ADJOINT:
3441 HU_inc = HU_inc + dESdLQ*ES_inc/HU_trl
3442 P_inc = P_inc + dESdP *ES_inc
3443 TT_inc = TT_inc + dESdTT*ES_inc
3444
3445 end subroutine HUtoES_ad
3446
3447 !--------------------------------------------------------------------------
3448 ! oop_calcGPSROJacobian
3449 !--------------------------------------------------------------------------
3450 subroutine oop_calcGPSROJacobian( columnTrlOnAnlIncLev, obsSpaceData, initializeLinearization_opt )
3451 !
3452 !:Purpose: Calculating the Jacobians of refractivity for oop_Hro/oop_HTro
3453 implicit none
3454
3455 ! Arguments:
3456 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
3457 type(struct_obs) , intent(inout) :: obsSpaceData
3458 logical, optional, intent(in) :: initializeLinearization_opt
3459
3460 ! Locals:
3461 real(8) :: zlat, lat
3462 real(8) :: zlon, lon
3463 real(8) :: zazm, azm
3464 integer :: isat
3465 real(8) :: rad, geo, zp0
3466 real(8), allocatable :: zpp(:), ztt(:), zhu(:), zHeight(:), zuu(:), zvv(:)
3467 real(8) :: zmt
3468 integer :: IDATYP, varNum
3469 integer :: jl, ngpslev, nwndlev
3470 integer :: headerIndex, bodyIndex, iProfile
3471 logical :: ASSIM
3472 logical, save :: initializeLinearization = .true.
3473 integer :: nh, nh1
3474 type(gps_profile) :: prf
3475 real(8) , allocatable :: h (:),azmv(:)
3476 type(gps_diff), allocatable :: rstv(:)
3477
3478 ! Re-compute the Jacobian for re-linearized state
3479 if ( present(initializeLinearization_opt) ) then
3480 initializeLinearization = initializeLinearization_opt
3481 end if
3482
3483 if ( .not. initializeLinearization ) return
3484
3485 write(*,*) 'ENTER oop_calcGPSROJacobian'
3486
3487 initializeLinearization = .false.
3488
3489 ! Initializations
3490 ngpslev = col_getNumLev(columnTrlOnAnlIncLev,'TH')
3491 nwndlev = col_getNumLev(columnTrlOnAnlIncLev,'MM')
3492
3493 allocate(zpp (ngpslev))
3494 allocate(ztt (ngpslev))
3495 allocate(zhu (ngpslev))
3496 allocate(zHeight (ngpslev))
3497 allocate(zuu (ngpslev))
3498 allocate(zvv (ngpslev))
3499
3500 if ( .not. allocated(oop_vRO_Jacobian) ) then
3501 allocate( oop_vRO_Jacobian(gps_numroprofiles,gps_ro_maxprfsize,4*ngpslev) )
3502 allocate( oop_vRO_lJac (gps_numROProfiles) )
3503 oop_vRO_Jacobian = 0.d0
3504 oop_vRO_lJac = .False.
3505 end if
3506
3507 allocate( h (gps_ro_maxprfsize) )
3508 allocate( azmv (gps_ro_maxprfsize) )
3509 allocate( rstv (gps_ro_maxprfsize) )
3510
3511 ! Loop over all header indices of the 'RO' family (Radio Occultation)
3512 ! Set the header list (start at the beginning of the list)
3513 call obs_set_current_header_list(obsSpaceData,'RO')
3514 HEADER: do
3515 headerIndex = obs_getHeaderIndex(obsSpaceData)
3516 if (headerIndex < 0) exit HEADER
3517
3518 ! Process only refractivity data (codtyp 169)
3519 IDATYP = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
3520 DATYP: if ( IDATYP == 169 ) then
3521 ! Scan for requested data values of the profile, and count them
3522 assim = .false.
3523 nh = 0
3524
3525 ! Loop over all body indices for this headerIndex:
3526 ! (start at the beginning of the list)
3527 call obs_set_current_body_list(obsSpaceData, headerIndex)
3528 BODY: do
3529 bodyIndex = obs_getBodyIndex(obsSpaceData)
3530 if (bodyIndex < 0) exit BODY
3531 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
3532 assim = .true.
3533 nh = nh + 1
3534 endif
3535 enddo BODY
3536
3537 ! If assimilations are requested, prepare and apply the observation operator
3538 ASSIMILATE: if (assim) then
3539 iProfile = gps_iprofile_from_index(headerIndex)
3540 if (oop_vRO_lJac(iProfile)) cycle ! If already done, end this HEADER
3541 varNum = gps_vRO_IndexPrf(iProfile, 2)
3542
3543 ! Profile at the observation location:
3544 ! Basic geometric variables of the profile:
3545 zlat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
3546 zlon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
3547 isat = obs_headElem_i(obsSpaceData,OBS_SAT,headerIndex)
3548 rad = obs_headElem_r(obsSpaceData,OBS_TRAD,headerIndex)
3549 geo = obs_headElem_r(obsSpaceData,OBS_GEOI,headerIndex)
3550 azm = obs_headElem_r(obsSpaceData,OBS_AZA,headerIndex)
3551 zmt = col_getHeight(columnTrlOnAnlIncLev,0,headerIndex,'SF')
3552 lat = zlat * MPC_DEGREES_PER_RADIAN_R8
3553 lon = zlon * MPC_DEGREES_PER_RADIAN_R8
3554 zazm = azm / MPC_DEGREES_PER_RADIAN_R8
3555 zp0 = col_getElem(columnTrlOnAnlIncLev,1,headerIndex,'P0')
3556 do jl = 1, ngpslev
3557 ! Profile x_b
3558 zpp(jl) = col_getPressure(columnTrlOnAnlIncLev,jl,headerIndex,'TH')
3559 ztt(jl) = col_getElem(columnTrlOnAnlIncLev,jl,headerIndex,'TT') - MPC_K_C_DEGREE_OFFSET_R8
3560 zhu(jl) = col_getElem(columnTrlOnAnlIncLev,jl,headerIndex,'HU')
3561 zHeight(jl) = col_getHeight(columnTrlOnAnlIncLev,jl,headerIndex,'TH')
3562 enddo
3563
3564 if ((col_getPressure(columnTrlOnAnlIncLev,1,headerIndex,'TH') + 1.0d-4) < &
3565 col_getPressure(columnTrlOnAnlIncLev,1,headerIndex,'MM')) then
3566 ! case with top thermo level above top momentum level (Vcode=5002)
3567 do jl = 1, nwndlev
3568 zuu(jl) = col_getElem(columnTrlOnAnlIncLev,jl,headerIndex,'UU')
3569 zvv(jl) = col_getElem(columnTrlOnAnlIncLev,jl,headerIndex,'VV')
3570 enddo
3571 else
3572 ! case without top thermo above top momentum level or unstaggered (Vcode=5001/4/5)
3573 do jl = 1, nwndlev-1
3574 zuu(jl) = col_getElem(columnTrlOnAnlIncLev,jl+1,headerIndex,'UU')
3575 zvv(jl) = col_getElem(columnTrlOnAnlIncLev,jl+1,headerIndex,'VV')
3576 enddo
3577 zuu(nwndlev) = zuu(nwndlev-1)
3578 zvv(nwndlev) = zuu(nwndlev-1)
3579 endif
3580 zuu(ngpslev) = zuu(nwndlev)
3581 zvv(ngpslev) = zuu(nwndlev)
3582
3583 ! GPS profile structure:
3584 call gps_struct1sw_v2(ngpslev,zlat,zlon,zazm,zmt,rad,geo,zp0,zpp,ztt,zhu,zuu,zvv,zHeight,prf)
3585
3586 ! Prepare the vector of all the observations:
3587 nh1 = 0
3588 call obs_set_current_body_list(obsSpaceData, headerIndex)
3589 BODY_2: do
3590 bodyIndex = obs_getBodyIndex(obsSpaceData)
3591 if (bodyIndex < 0) exit BODY_2
3592 if ( obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
3593 nh1 = nh1 + 1
3594 h(nh1) = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
3595 azmv(nh1)= zazm
3596 endif
3597 enddo BODY_2
3598
3599 ! Apply the observation operator:
3600 ! varNum = bufr_nebd (15037) or varNum = bufr_nerf (15036) for GPS-RO
3601 if (varNum == bufr_nebd) then
3602 call gps_bndopv1(h, azmv, nh, prf, rstv)
3603 else
3604 call gps_refopv (h, nh, prf, rstv)
3605 end if
3606 do nh1 = 1, nh
3607 oop_vRO_Jacobian(iprofile,nh1,1:4*ngpslev)= rstv(nh1)%dvar(1:4*ngpslev)
3608 end do
3609 oop_vRO_lJac(iProfile) = .True.
3610 endif ASSIMILATE
3611 endif DATYP
3612 enddo HEADER
3613
3614 deallocate( rstv )
3615 deallocate( azmv )
3616 deallocate( h )
3617
3618 deallocate(zvv)
3619 deallocate(zuu)
3620 deallocate(zHeight)
3621 deallocate(zhu)
3622 deallocate(ztt)
3623 deallocate(zpp)
3624
3625 write(*,*) 'EXIT oop_calcGPSROJacobian'
3626
3627
3628 end subroutine oop_calcGPSROJacobian
3629
3630 !--------------------------------------------------------------------------
3631 ! oop_calcGPSGBJacobian
3632 !--------------------------------------------------------------------------
3633 subroutine oop_calcGPSGBJacobian( columnTrlOnAnlIncLev, obsSpaceData, initializeLinearization_opt )
3634 !
3635 !:Purpose: Calculating the Jacobians of ZTD for oop_Hgp/oop_HTgp
3636 implicit none
3637
3638 ! Arguments:
3639 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
3640 type(struct_obs) , intent(inout) :: obsSpaceData
3641 logical, optional , intent(in) :: initializeLinearization_opt
3642
3643 ! Locals:
3644 real(8) :: ZLAT, Lat
3645 real(8) :: ZLON, Lon
3646 real(8), allocatable :: ZTTB(:)
3647 real(8), allocatable :: ZHUB(:)
3648 real(8), allocatable :: zHeight(:)
3649 real(8), allocatable :: ZPPB(:)
3650 real(8) :: ZP0B, ZPSMOD, ZPWMOD, ZPWMOD2, dZTD
3651 real(8) :: ZMT
3652 real(8) :: sfcfield
3653 real(8) :: dxq1, dxq2, dxq3
3654 real(8) :: ZLEV, ZDZMIN
3655 real(8) :: JAC(gps_ncvmx)
3656 integer :: headerIndex, bodyIndex
3657 integer :: JL, NFLEV, iztd, icount, vcode
3658 logical :: ASSIM
3659 type(gps_profilezd) :: PRF, PRF2
3660 type(gps_diff) :: ZTDOPV, ZTDOPV2
3661 type(struct_vco), pointer :: vco_anl
3662 character(len=12) :: cstnid
3663 logical, save :: initializeLinearization = .true.
3664
3665 ! Re-compute the Jacobian for re-linearized state
3666 if ( present(initializeLinearization_opt) ) then
3667 initializeLinearization = initializeLinearization_opt
3668 end if
3669
3670 if ( .not. initializeLinearization ) return
3671
3672 write(*,*) 'ENTER oop_calcGPSGBJacobian'
3673
3674 initializeLinearization = .FALSE.
3675
3676 vco_anl => col_getVco(columnTrlOnAnlIncLev)
3677 vcode = vco_anl%vCode
3678
3679 ZDZMIN = gps_gb_DZMIN ! from modgpsztd_mod
3680
3681 NFLEV = col_getNumLev(columnTrlOnAnlIncLev,'TH')
3682
3683 ! Initializations
3684 if ( .not. allocated(gps_ZTD_Index) ) call utl_abort('oop_calcGPSGBJacobian: ERROR: gps_ZTD_Index not allocated!')
3685 if ( allocated(oop_vZTD_Jacobian) ) write(*,*) 'oop_calcGPSGBJacobian: WARNING, oop_vZTD_Jacobian is already allocated. Re-allocating'
3686 call utl_reallocate(oop_vZTD_Jacobian,gps_gb_numZTD,4*NFLEV)
3687 oop_vZTD_Jacobian(:,:) = 0.0d0
3688
3689 allocate(ZTTB(NFLEV))
3690 allocate(ZHUB(NFLEV))
3691 allocate(zHeight(NFLEV))
3692 allocate(ZPPB(NFLEV))
3693
3694 write(*,*) 'oop_calcGPSGBJacobian: Storing Jacobians for GPS ZTD data ...'
3695 write(*,*) ' INFO: Analysis grid iversion = ', vcode
3696 write(*,*) ' col_getNumLev[columnTrlOnAnlIncLev,TH] = ', NFLEV
3697 write(*,*) ' gps_gb_numZTD = ', gps_gb_numZTD
3698
3699 icount = 0
3700
3701 ! loop over all header indices of the 'GP' family (GPS observations)
3702 call obs_set_current_header_list(obsSpaceData,'GP')
3703 HEADER_0: do
3704 headerIndex = obs_getHeaderIndex(obsSpaceData)
3705 if (headerIndex < 0) exit HEADER_0
3706
3707 ! Scan for ZTD assimilation at this location
3708 ASSIM = .FALSE.
3709
3710 ! loop over all body indices for this headerIndex
3711 call obs_set_current_body_list(obsSpaceData, headerIndex)
3712 BODY_0: do
3713 bodyIndex = obs_getBodyIndex(obsSpaceData)
3714 if (bodyIndex < 0) exit BODY_0
3715 if ( (obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex) == 189) &
3716 .and. (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == bufr_nezd) &
3717 .and. (obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) ) then
3718 ZLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
3719 ASSIM = .TRUE.
3720 end if
3721 end do BODY_0
3722
3723 if ( ASSIM ) THEN
3724 ! LR background profile at the observation location x :
3725 icount = icount + 1
3726 Lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
3727 ZLAT = Lat * MPC_DEGREES_PER_RADIAN_R8
3728 Lon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
3729 ZLON = Lon * MPC_DEGREES_PER_RADIAN_R8
3730 ZP0B = col_getElem(columnTrlOnAnlIncLev,1,headerIndex,'P0')
3731 do JL = 1, NFLEV
3732 ZTTB(JL) = col_getElem(columnTrlOnAnlIncLev,JL,headerIndex,'TT') - MPC_K_C_DEGREE_OFFSET_R8
3733 ZHUB(JL) = col_getElem(columnTrlOnAnlIncLev,JL,headerIndex,'HU')
3734 ZPPB(JL) = col_getPressure(columnTrlOnAnlIncLev,JL,headerIndex,'TH')
3735 zHeight(JL) = col_getHeight(columnTrlOnAnlIncLev,JL,headerIndex,'TH')
3736 end do
3737 if ( abs(ZPPB(NFLEV)-ZP0B) > 0.1 ) then
3738 write(*,*) ' oop_calcGPSGBJacobian: ERROR: |ZPPB(NFLEV)-ZP0B| > 0.1'
3739 write(*,*) ' ZPPB(NFLEV), ZP0B =', ZPPB(NFLEV), ZP0B
3740 end if
3741 ZMT = col_getHeight(columnTrlOnAnlIncLev,0,headerIndex,'SF')
3742
3743 call gps_structztd_v2(NFLEV,Lat,Lon,ZMT,ZP0B,ZPPB,ZTTB,ZHUB,zHeight,gps_gb_LBEVIS,gps_gb_IREFOPT,PRF)
3744 call gps_ztdopv(ZLEV,PRF,gps_gb_LBEVIS,ZDZMIN,ZTDopv,ZPSMOD,gps_gb_IZTDOP)
3745
3746 ! Observation Jacobian H'(xb)
3747 JAC = ZTDopv%DVar
3748 iztd = gps_iztd_from_index(headerIndex)
3749 do JL = 1, 4*NFLEV
3750 oop_vZTD_Jacobian(iztd,JL) = JAC(JL)
3751 end do
3752
3753 if ( icount <= 3 .and. gps_gb_LTESTOP ) then
3754 write(*,*) '--------------------------------------------------------- '
3755 cstnid = obs_elem_c(obsSpaceData,'STID',headerIndex)
3756 write(*,*) iztd, cstnid, 'ZTDopv (m) = ', ZTDopv%Var
3757 call gps_pw(PRF,ZPWMOD)
3758
3759 sfcfield = ZP0B + 50.0d0
3760 call gps_structztd_v2(NFLEV,Lat,Lon,ZMT,sfcfield,ZPPB,ZTTB,ZHUB,zHeight,gps_gb_LBEVIS,gps_gb_IREFOPT,PRF2)
3761 call gps_ztdopv(ZLEV,PRF2,gps_gb_LBEVIS,ZDZMIN,ZTDopv2,ZPSMOD,gps_gb_IZTDOP)
3762 write(*,*) ' ZTD Operator Test: dP0 = +50 Pa'
3763 write(*,*) ' dZTD NL = ', ZTDopv2%Var - ZTDopv%Var
3764 write(*,*) ' dZTD Linear = ', oop_vZTD_Jacobian(iztd,4*NFLEV)*50.0d0
3765 write(*,*) ' '
3766
3767 ! q dx
3768 dxq1 = 0.44D-01*ZHUB(64)
3769 dxq2 = 0.44D-01*ZHUB(65)
3770 dxq3 = 0.44D-01*ZHUB(66)
3771 ZHUB(64) = ZHUB(64) - dxq1
3772 ZHUB(65) = ZHUB(65) - dxq2
3773 ZHUB(66) = ZHUB(66) - dxq3
3774 call gps_structztd_v2(NFLEV,Lat,Lon,ZMT,ZP0B,ZPPB,ZTTB,ZHUB,zHeight,gps_gb_LBEVIS,gps_gb_IREFOPT,PRF2)
3775 call gps_ztdopv(ZLEV,PRF2,gps_gb_LBEVIS,ZDZMIN,ZTDopv2,ZPSMOD,gps_gb_IZTDOP)
3776 call gps_pw(PRF2,ZPWMOD2)
3777 write(*,*) ' ZTD Operator Test: dQ = -0.44E-01*Q JL = 64,65,66'
3778 write(*,*) ' dPW (mm) = ', ZPWMOD2 - ZPWMOD
3779 write(*,*) ' dZTD NL = ', ZTDopv2%Var - ZTDopv%Var
3780 dZTD = oop_vZTD_Jacobian(iztd,64+NFLEV)*(-dxq1) + oop_vZTD_Jacobian(iztd,65+NFLEV)*(-dxq2) + &
3781 oop_vZTD_Jacobian(iztd,66+NFLEV)*(-dxq3)
3782 write(*,*) ' dZTD Linear = ', dZTD
3783 write(*,*) ' '
3784 ZHUB(64) = ZHUB(64) + dxq1
3785 ZHUB(65) = ZHUB(65) + dxq2
3786 ZHUB(66) = ZHUB(66) + dxq3
3787
3788 ! temperature dx
3789 ZTTB(64) = ZTTB(64) + 2.0d0
3790 ZTTB(65) = ZTTB(65) + 2.0d0
3791 ZTTB(66) = ZTTB(66) + 2.0d0
3792 call gps_structztd_v2(NFLEV,Lat,Lon,ZMT,ZP0B,ZPPB,ZTTB,ZHUB,zHeight,gps_gb_LBEVIS,gps_gb_IREFOPT,PRF2)
3793 call gps_ztdopv(ZLEV,PRF2,gps_gb_LBEVIS,ZDZMIN,ZTDopv2,ZPSMOD,gps_gb_IZTDOP)
3794 write(*,*) ' ZTD Operator Test: dTT = +2.0K JL = 64,65,66'
3795 write(*,*) ' dZTD NL = ', ZTDopv2%Var - ZTDopv%Var
3796 dZTD = oop_vZTD_Jacobian(iztd,64)*2.0d0 + oop_vZTD_Jacobian(iztd,65)*2.0d0 + &
3797 oop_vZTD_Jacobian(iztd,66)*2.0d0
3798 write(*,*) ' dZTD Linear = ', dZTD
3799 write(*,*) '--------------------------------------------------------- '
3800 end if
3801
3802 end if
3803
3804 end do HEADER_0
3805
3806 deallocate(ZTTB)
3807 deallocate(ZHUB)
3808 deallocate(zHeight)
3809 deallocate(ZPPB)
3810
3811 write(*,*) 'oop_calcGPSGBJacobian: Number of ZTD data (icount) = ', icount
3812 write(*,*) ' Expected number (gps_gb_numZTD) = ', gps_gb_numZTD
3813 write(*,*) ' Last iztd = ', iztd
3814 write(*,*) ' gps_ZTD_Index(1) = ', gps_ZTD_Index(1)
3815 write(*,*) ' gps_ZTD_Index(iztd) = ', gps_ZTD_Index(iztd)
3816
3817 if ( icount /= gps_gb_numZTD ) then
3818 call utl_abort('oop_calcGPSGBJacobian: ERROR: icount /= gps_gb_numZTD!')
3819 end if
3820 if ( icount /= iztd ) then
3821 call utl_abort('oop_calcGPSGBJacobian: ERROR: icount /= iztd!')
3822 end if
3823 if ( gps_gb_numZTD /= iztd ) then
3824 call utl_abort('oop_calcGPSGBJacobian: ERROR: gps_gb_numZTD /= iztd!')
3825 end if
3826
3827 write(*,*) 'EXIT oop_calcGPSGBJacobian'
3828
3829 end subroutine oop_calcGPSGBJacobian
3830
3831 function oop_iceScaling(obsSpaceData, bodyIndex) result(scaling)
3832 !
3833 !:Purpose: Calculate the scaling factor for
3834 ! ice related observations to convert from
3835 ! model space to observation space, i.e.
3836 ! H(iceConc) = scaling*iceConc + constant
3837 !
3838 implicit none
3839
3840 ! Arguments:
3841 type(struct_obs), intent(in) :: obsSpaceData
3842 integer, intent(in) :: bodyIndex
3843 ! Result:
3844 real(8) :: scaling
3845
3846 ! Locals:
3847 integer :: bufrCode
3848 integer :: headerIndex, obsDate, monthIndex
3849 integer :: trackCellNum
3850 character(len=8) :: ccyymmdd
3851
3852 bufrCode = obs_bodyElem_i( obsSpaceData, OBS_VNM, bodyIndex )
3853
3854 select case (bufrCode)
3855 case(BUFR_ICEC, BUFR_ICEP)
3856 scaling = 100.0d0
3857 case(BUFR_ICEV)
3858 scaling = 1.0d0
3859 case(BUFR_ICES)
3860 headerIndex = obs_bodyElem_i( obsSpaceData, OBS_HIND, bodyIndex )
3861 obsDate = obs_headElem_i( obsSpaceData, OBS_DAT, headerIndex )
3862 write(ccyymmdd, FMT='(i8.8)') obsDate
3863 read(ccyymmdd(5:6), FMT='(i2)') monthIndex
3864 trackCellNum = obs_headElem_i( obsSpaceData, OBS_FOV, headerIndex )
3865 scaling = oer_ascatAnisIce(trackCellNum,monthIndex) - &
3866 oer_ascatAnisOpenWater(trackCellNum,monthIndex)
3867 case default
3868 scaling = 0.0d0
3869 end select
3870
3871 end function oop_iceScaling
3872
3873end module obsOperators_mod
