1module obsOperatorsChem_mod
2 ! MODULE obsOperatorsChem_mod (prefix='oopc' category='5. Observation operators')
3 !
4 !:Purpose: Observation operators for CH obs family, including nonlinear, tangent-linear
5 ! and adjoint versions, and related setup and input routines.
6 !
7 use earthConstants_mod
8 use mathPhysConstants_mod
9 use obsSpaceData_mod
10 use columnData_mod
11 use bufr_mod
12 use physicsFunctions_mod
13 use midasMpi_mod
14 use utilities_mod
15 use varNameList_mod
16 use obsSubSpaceData_mod
17 use obsFiles_mod
18 use codtyp_mod
19 use ozoneClim_mod
20 use presProfileOperators_mod
21 use timeCoord_mod
22 use bCovarSetupChem_mod
23
24 implicit none
25 save
26 private
27
28 ! public procedures
29 public :: oopc_CHobsoperators, oopc_diagnOnly, oopc_addEfftempObsfile
30
31 !-------------------------------------------------------------------------
32 ! Various structures and parameters for the CH family
33
34 ! module structures
35 ! -----------------
36
37 type :: struct_oopc_obsoperators
38
39 ! Structure holding work variables for individual observation operator
40 ! applications
41 !
42 ! Variable Description
43 ! -------- -----------
44 ! nobslev Number of observations in the profile
45 ! nmodlev Number of model levels in the column
46 ! varno BUFR descriptor element for obs units
47 ! constituentId BUFR code element of local GRIB Table 08046 identifying the constituent
48 ! (similar to BUFR Table 08043)
49 ! operatorCategory CH obs operator category. One of the following.
50 ! 'Interp': interpolators
51 ! 'Surface': surface point operators
52 ! 'Integ': integration operators
53 ! 'LayerAvg': layer averaging operators
54 ! Automatically assigned based on obs BUFR element
55 ! and, when relevant, assocciated obsinfo_chm content.
56 ! layerIdentified .true. if a layer (with identified layer boundaries)
57 ! .false if layer boundaries are not available.
58 ! valertop Layer top (final work values in Pa)
59 ! vlayerbottom Layer bottom (final work values in Pa)
60 ! zh Initial innovation model array (other than conversion constants)
61 ! zhp Part of innovation operator not related to resolution.
62 ! modlevindexTop Top level of non-zero values in zh
63 ! modlevindexBot Bottom level of non-zero values in zh
64 ! trial Trial (background) profile on model levels at observation location
65 ! tt Temperature profile on model levels (Kelvin)
66 ! hu Specific humidity on model levels
67 ! height Height on model levels (m)
68 ! pp Pressure on model levels (Pa)
69 ! lat Latitude of observation (radians)
70 ! lon Longitude of observation (radians)
71 ! obslev Observation profile level values (OBS_PPP)
72 ! obsSpaceTrial obs estimate on obs levels based on trial/background profile
73 ! varName Variable/obs nomvar
74 ! stnid Observation station ID
75 ! date YYYYMMDD (date of obs)
76 ! hhmm HHMM (time of obs)
77 ! obs_index Observation index
78 ! Note: Depending on the data of interest, the index of a required array element or
79 ! profile associated to an observation can be identified from (lat,long,date,hhmm,
80 ! stnid, optional task-dependent identifier if needed) or "obs_index".
81 ! The latter is for associations of data identified within
82 ! processing of individual CPUs. Each of the two index identifiers is represented
83 ! by the unique character string identifier 'code' of struct_oss_obsdata
84 ! (e.g. see obsdata_get_header_code_r for use of (lat,long,date,hhmm,stnid)).
85 ! vco Index of vertical coord type for obs
86 ! 1 - Altitudes (m)
87 ! 2 - Pressure (Pa)
88 ! 3 - Channel index
89 ! 4 - not provided with obs. Obs is for total column values.
90 ! 5 - not provided with obs. Obs is a surface point value.
91 ! iavgkern Integer indicating if averaging kernels are to be applied. Value
92 ! of zero indicates no averaging kernel to be applied. Non-zero value
93 ! indicates index in oopc_avgkern%obsSubSpace arrays.
94 ! applyGenOper Indicates if the generalized observation operator should be applied
95 ! columnBound Boundary imposed on a column measurement
96 ! ixtr Indicates when outside model vertical range (when >0)
97 ! success Indicates if the observation is to be assimilated
98 ! or was successfully assimilated
99
100 integer :: nobslev,nmodlev,constituentId,vco,varno,date,hhmm,iavgkern,obs_index
101 logical :: layerIdentified,applyGenOper
102 real(8) :: lat,lon,columnBound
103 character(len=12) :: stnid
104 character(len=10) :: operatorCategory
105 character(len=4) :: varName
106 real(8), allocatable :: vlayertop(:),vlayerbottom(:),tt(:),height(:),pp(:)
107 real(8), allocatable :: zh(:,:),zhp(:,:),obslev(:),hu(:)
108 integer, allocatable :: ixtr(:)
109 logical, allocatable :: success(:)
110 real(8), pointer :: trial(:)
111 integer, allocatable :: modlevindexTop(:),modlevindexBot(:)
112 real(8), allocatable :: obsSpaceTrial(:)
113
114 end type struct_oopc_obsoperators
115
116 type(struct_oopc_obsoperators) :: obsoper
117
118 type :: struct_oopc_operatorsDepot
119
120 ! Structure holding saved arrays for observation operators
121 ! A subset of 'struct_oopc_obsoperators'
122 !
123 ! Variable Description
124 ! -------- -----------
125 ! nobslev Number of observations in the profile
126 ! valertop Layer top (final work values in Pa)
127 ! vlayerbottom Layer bottom (final work values in Pa)
128 ! zh Initial innovation model array (other than conversion constants)
129 ! zhp Part of innovation operator not related to resolution.
130 ! modlevindexTop Top level of non-zero values in zh
131 ! modlevindexBot Bottom level of non-zero values in zh
132 ! trial Trial (background) profile on model levels at observation location
133 ! iavgkern Integer indicating if averaging kernels are to be applied. Value
134 ! of zero indicates no averaging kernel to be applied. Non-zero value
135 ! indicates index in oopc_avgkern%obsSubSpace arrays.
136 ! applyGenOper Indicates if the generalized observation operator should be applied
137 ! ixtr Indicates when outside model vertical range (when >0)
138 ! success Indicates if the observation is to be assimilated
139 ! or was successfully assimilated
140
141 integer :: nobslev,iavgkern
142 logical :: applyGenOper
143 real(8), allocatable :: vlayertop(:),vlayerbottom(:)
144 real(8), allocatable :: zh(:,:),zhp(:,:)
145 integer, allocatable :: ixtr(:)
146 logical, allocatable :: success(:)
147 real(8), pointer :: trial(:)
148 integer, allocatable :: modlevindexTop(:),modlevindexBot(:)
149
150 end type struct_oopc_operatorsDepot
151
152 type :: struct_oopc_info
153
154 ! Information arrays retrieved from auxiliary file regarding vertical levels
155 ! or averaging kernels
156 !
157 ! Variable Description
158 ! -------- -----------
159 ! n_stnid Number of sub-families (identified via STNIDs)
160 ! stnids Sub-families (STNIDs; * are wild cards)
161 ! element BUFR element in data block
162 ! profElement BUFR element of profile required by obs operator if
163 ! if it differs from 'element' (e.g. see 'n_col' and 'type' below).
164 ! source 0: Set entirely from the auxiliary file being read. No
165 ! initial values read from observation files
166 ! 1: Initial values in observation files for constant number
167 ! of vertical levels (may be adjusted after input)
168 ! vco Index of vertical coord type for obs
169 ! 1 - Altitudes (m)
170 ! 2 - Pressure (Pa)
171 ! 3 - Channel index
172 ! 4 - not provided with obs. Obs is for total column values.
173 ! 5 - not provided with obs. Obs is a surface point value.
174 ! 6 - sigma vertical coordinate (*Psfc to give levels in Pa); in this module only.
175 ! ibegin Position index of start of data for given
176 ! sub-family.
177 ! n_lvl Lengh of corresponding obs profile (number of rows)
178 ! n_col Length of array for each obs element (number of columns; optional)
179 ! - Usually n_lvl will be the number of vertical levels of profiles
180 ! and n_col is not needed.
181 ! - For averaging kernels:
182 ! - When the obs are single level data, n_lvl=1 or 2. If it requires more levels
183 ! or data for observation operators, n_col = number of additional data
184 ! - The a-priori contribution (I-A)xa can be added as the element on each row.
185 ! When done and n_lvl>1, n_col would be expected to be n_lvl+1
186 ! - If an integeration over the n_lvl values (>1) needs to be performed,
187 ! the column n_col = n_lvl+2 provides the integration weights for a summations
188 ! that gives the integral.
189 ! - For vertical levels:
190 ! - n_col will be set to 2 as default for layer boundaries
191 ! - n_col needs to be set to 1 to use this simply as a single vertical level profile
192 ! n_lat Number of latitudes
193 ! lat Latitudes (degrees; ordered in increasing size)
194 !
195 ! vlayertop Layer top
196 ! vlayerbottom Layer bottom
197 ! rak Averaging kernel matrices
198 ! type Operation subtypes.
199 !
200 ! Currently relevant only for averaging kernels:
201 !
202 ! Type name Description
203 ! ============================================================
204 ! 'default' No special treatment (default)
205 ! 'log' Application of log-space averaging kernel
206 ! ============================================================
207 !
208 integer :: n_stnid
209 character(len=12), allocatable :: stnids(:)
210 character(len=15), allocatable :: type(:)
211 integer, allocatable :: element(:),source(:),profElement(:)
212 integer, allocatable :: vco(:),n_lat(:)
213 integer, allocatable :: ibegin(:),n_lvl(:),n_col(:)
214 real(8), allocatable :: rak(:),vlayertop(:),vlayerbottom(:)
215 real(8), allocatable :: lat(:)
216
217 type(struct_oss_obsdata), allocatable :: obsSubSpace(:)
218
219 end type struct_oopc_info
220
221 type(struct_oopc_info) :: oopc_levels
222 type(struct_oopc_info) :: oopc_avgkern
223
224 ! Arrays for integration upper boundary of retrieved total column measurements
225 type(struct_oss_obsdata) :: oopc_columnBoundary
226
227 ! File name of auxiliary text file constaining supplemental observation information
228 character(len=50), parameter :: oopc_aux_filename="obsinfo_chm"
229
230 ! Max nummber of constituents (max size of related arrays)
231 integer, parameter :: oopc_constituentsSize=30 ! = max allowed value of "iconstituentId" for Table 08046.
232 ! Value to be increased as needed up to a max of 6999 as values
233 ! > 7000 (and less 0) are assumed assigned to non-constituent fields
234
235 ! Arrays containing input reference fields and fields interpolated
236 ! to obs locations
237
238 type :: struct_oopc_field
239
240 ! Structure for storing reference (climatological) fields needed for
241 ! operatorSubType(2,:) == 'genOper' with
242 ! genOperConstraintType == 'Diff' (see below).
243 !
244 ! Variable Description
245 ! -------- -----------
246 ! field Gridded 3D field (lon,lat,vlev) or 2D field (1,lat,vlev)
247 ! nlat number of latitudes
248 ! nlon number of longitudes
249 ! nlev number of vertical levels
250 ! lat,lon lat,lon grid in radians
251 ! vlev vertical levels
252 ! ivkind Index of vertical coordinate type. Defintion may vary according to source.
253 ! For fields read from RPN files and use of convip:
254 ! 0: P is in height [m] (metres) with respect to sea level
255 ! 1: P is in stddev [sg] (0.0 -> 1.0)
256 ! 2: P is in pressure [mb] (millibars)
257 ! 3: P is in an arbitrary code
258 ! 4: P is in height [M] (metres) with respect to ground level
259 ! 5: P is in hybrid coordinates [hy]
260 ! 6: P is in theta [th]
261 ! For use with obs
262
263 real(8), allocatable :: field(:,:,:),lat(:),lon(:),vlev(:)
264 integer :: nlev,nlon,nlat,ivkind
265
266 end type struct_oopc_field
267
268 type(struct_oss_obsdata) :: oopc_bgRef
269 type(struct_oopc_field) :: oopc_climFields(0:oopc_constituentsSize,2)
270
271 ! Arrays to contain the calculated concentration-weighted effective temperature
272 ! associated to total column data. It will be stored in the observation file.
273
274 type(struct_oss_obsdata) :: oopc_efftemp
275
276 type(struct_bcsc_bgStats) :: bgStats ! Background covariances
277
278 ! Setup initialization key
279 logical :: initializedChem = .false.
280
281 ! Operator storage initialization key
282 logical :: initializedOperators = .false.
283
284 ! General config/setup information parameters
285 ! See description list of NAMCHEM namelist parameters and others
286 ! in routine oopc_readNamchem.
287 ! Following variables/parameters could be placed in a data structure/type
288 ! (e.g. struct_oopc_nmlparm)
289
290 character(len=5) :: oopc_genOperConstraintType(0:oopc_constituentsSize)
291 real(8) :: oopc_genOperHCorrlenExpnt(0:oopc_constituentsSize)
292 real(8) :: oopc_genOperOmAStatsFactor(0:oopc_constituentsSize)
293 integer :: oopc_tropo_mode(0:oopc_constituentsSize),oopc_tropo_bound(0:oopc_constituentsSize)
294 integer :: oopc_obsdata_maxsize
295 real(8) :: oopc_tropo_column_top(0:oopc_constituentsSize)
296 logical :: oopc_storeOperators
297
298 integer, parameter :: assim_maxfamnum=1 ! Could be used for other families as well with >1
299 integer, parameter :: assim_maxsize=100 ! max size of assim_* arrays
300
301 integer :: assim_famNum
302 character(len=2) :: assim_fam(assim_maxfamnum) ! List of families to which filt_diagnOnly is to apply
303 logical :: assim_all(assim_maxfamnum) ! Choose to assimilate all obs of specified family
304 integer :: assim_num(assim_maxfamnum) ! Number of combinations identified for assimilation
305 integer :: assim_varno(assim_maxfamnum,assim_maxsize) ! List of bufr elements to assimilate (0 means all)
306 integer :: assim_nlev(assim_maxfamnum,assim_maxsize) ! 0: multi- and uni-lev; 1: uni-lev; >1 multi-lev
307 integer :: assim_exclude_nflag(assim_maxfamnum) ! List of bits for excluding obs from assimilation
308 integer :: assim_exclude_flag(assim_maxfamnum,assim_maxsize) ! Number of bits for excluding obs
309 character(len=9) :: assim_stnid(assim_maxfamnum,assim_maxsize) ! List of stnid to assimilation '*' for wild card
310 character(len=20) :: operatorSubType(2,assim_maxsize) ! Operator sub-type name
311 character(len=10) :: modelName = 'GEM-MACH' ! Identification of the model
312
313 !**************************************************************************
314
315 contains
316
317 !--------------------------------------------------------------------------
318 ! oopc_setupCH
319 !--------------------------------------------------------------------------
320 subroutine oopc_setupCH(kmode)
321 !
322 !:Purpose: To set up additional information required by constituent obs and
323 ! not provided in obsSpaceData. Also to assign observation layer
324 ! top and bottom levels (and averaging kernel matrices).
325 ! See 'oopc_CHobsoperators'.
326 !
327 implicit none
328
329 ! Arguments:
330 integer, intent(in) :: kmode ! Mode of observation operator
331
332 ! Locals:
333 logical :: success
334
335 write(*,*) 'Begin oopc_setupCH'
336
337 ! Read NAMCHEM namelist and set related parameters
338
339 call oopc_readNamchem
340 write(*,*) 'oopc_setupCH: Completed oopc_readNamchem'
341
342 ! Read reference vertical levels (where needed) or
343 ! top and bottom layer boundaries of partial (or total) column meausurements
344
345 call oopc_readLevels
346 write(*,*) 'oopc_setupCH: Completed oopc_readLevels'
347
348 ! To deallocate space if required elsewhere, one should use
349 ! call oopc_deallocLevels
350
351 ! Read averaging kernel matrices
352
353 call oopc_readAvgkern
354 write(*,*) 'oopc_setupCH: Completed oopc_readAvgkern'
355
356 ! To deallocate space if required elsewhere, one should use
357 ! call oopc_deallocAvgkern
358
359 ! Read reference (e.g. climatological) fields
360
361 call oopc_readFields(oopc_climFields,oopc_aux_filename,'CH', &
362 oopc_constituentsSize,2,oopc_genOperConstraintType, &
363 success,filetype_opt='TXT')
364 if ( .not. success ) then
365 call utl_abort('oopc_setupCH: Failed in oopc_readFields')
366 end if
367 write(*,*) 'oopc_setupCH: Completed oopc_readFields'
368
369 ! Allocation of oopc_efftemp done in oopc_setupCH instead of obsdata_add_data1d
370 ! to ensure allocation is done for all processors, including those without associated data.
371 ! This is to ensure that rpn_comm_allgather will work in routine obsdata_MPIGather.
372
373 if (.not.associated(oopc_efftemp%data1d)) then
374 call oss_obsdata_alloc(oopc_efftemp,oopc_obsdata_maxsize,dim1=1)
375 oopc_efftemp%nrep=0
376 end if
377
378 ! Get background error std dev stats when they may be required
379
380 bgStats%initialized = .false.
381
382 if (kmode == 1 .or. any(operatorSubType(2,:) == 'genOper') ) then
383 call bcsc_getCovarCH(bgStats)
384 end if
385
386 write(*,*) 'Completed oopc_setupCH'
387
388 end subroutine oopc_setupCH
389
390 !--------------------------------------------------------------------------
391 ! oopc_readNamchem
392 !--------------------------------------------------------------------------
393 subroutine oopc_readNamchem
394 !:Purpose: Read and store miscellaneous flags and constants.
395 !
396 !:Comment: assim_* arrays could instead be made available to all families
397 ! by moving them to a different input namelist (and changing its
398 ! dimensions settings).
399 !
400 !
401 ! :genOperConstraintType:
402 ! Reference profile type for weighted integration
403 ! or layer averaging (generalized observation)
404 ! operator.
405 ! Relevant for operatorSubType(2,i)='genOper'.
406 ! ================================================
407 ! 'Trial' use trial field xb for mass weighted
408 ! increment distribution
409 ! 'Diff' use a combination of the difference of
410 ! an external reference xc and the trial
411 ! field xb, i.e. mass weighted increment
412 ! distribution as a(xc-xb) + b*xc where a
413 ! and b depend on the size of
414 ! sum[(xc-xb)/sig(xb)]^2 over the profile
415 ! ==============================================
416 !
417 ! :genOperHCorrlenExpnt: Used with operatorSubType(2,i) ='genOper'
418 ! Exponent for partially mitigating the effect of
419 ! the influence of neighbouring column amonunt obs
420 ! from background error correlations.
421 ! Emperically obtained exponent value.
422 ! Not optimal for all possible local horizontal data densities.
423 !
424 ! :genOperOmAStatsFactor: OmA RMS (or std dev) conservation factor for
425 ! operatorSubType(2,i) ='genOper'.
426 !
427 !
428 ! :assim_fam: List of families to which filt_diagnOnly is to
429 ! apply.
430 !
431 ! :assim_exclude_flag: Array specifying bits for identifying
432 ! diagnostic-only observations for observations
433 ! that would otherwise be assimilated according to
434 ! the other assim_* arrays
435 !
436 ! :assim_exclude_nflag: Number of bit flags to specify in
437 ! assim_exclude_flag array
438 !
439 ! :assim_all: Logical indicating if all assimilatable obs of
440 ! the specified family will be assimilated
441 ! (default is .true.)
442 ! When assim_all is .false., account for the setttings
443 ! of assim_num, assim_varno, assim_stnid, assim_nlev.
444 !
445 ! :assim_num: Relevant when assim_all = ,false.
446 ! Number combinations (stnid, bufr element,
447 ! multi/uni-level) identified for assimilation.
448 ! All others will not be assimilated. OmP and OmA
449 ! diagnostics and output will still be produced
450 ! for non-assimilated datasets.
451 !
452 ! === =======================================
453 ! 0 none are to be assimilated when assim_all
454 ! is .false. (default)
455 ! >0 sets of (stnid, bufr varno,
456 ! multi/uni-levels) to be assimilated
457 ! === =======================================
458 !
459 ! :assim_varno: Bufr elements of obs sets for assimilation. A
460 ! value of 0 implies that all are to be used.
461 !
462 ! :assim_stnid: Stnids of obs sets for assimilation. '*' denote
463 ! wild cards
464 !
465 ! :assim_nlev: === =========================
466 ! 0 multi-level and uni-level
467 ! 1 uni_level
468 ! >1 multi-level
469 ! === =========================
470 !
471 ! :tropo_mode: Integer indicating if special treatment is to be
472 ! given to the troposphere when assimilating total
473 ! column measurements. Values indicate
474 ! === =======================================
475 ! 0 No special treatment given (default)
476 ! 1 Values of the adjoint model above
477 ! obsoper%columnBound set to zero. If
478 ! specified, generalized innovation
479 ! operator only applied below
480 ! obsoper%columnBound in the tangent
481 ! linear model.
482 ! 2 Values of tangent linear model and
483 ! adjoint model above obsoper%columnBound
484 ! set to zero.
485 ! === =======================================
486 ! Array index refers to BUFR code element of Table
487 ! 08046 (iconstituentId) identifying the
488 ! constituent. Relevant for total column
489 ! measurements only.
490 !
491 ! :tropo_bound: Integer indicating which column top value to use
492 ! if tropo_mode is non-zero.
493 ! === =======================================
494 ! 0 Use fixed value of tropo_column_top
495 ! 1 Use model determination of tropopause
496 ! 2 Use model determination of PBL
497 ! === =======================================
498 ! Options 1 and 2 will default to the value set
499 ! in tropo_column_top if the model derived column
500 ! top could not be determined. Relevant for total
501 ! column measurements only.
502 !
503 ! :tropo_column_top: Default value to use for the column boundary
504 ! (in Pa). Array index refers to BUFR code element
505 ! of Table 08046 (iconstituentId) identifying
506 ! the constituent. Relevant for total column
507 ! measurements only.
508 !
509 ! :obsdata_maxsize: Max allowed size of work arrays (in terms of
510 ! number of obs) associated to ordered observation
511 ! indices
512 !
513 !
514 ! :modelName: Identifier of forecast model
515 ! Default: 'GEM-MACH'
516 ! Set to 'GEM' for varNames of 'O3L', 'CH4L', and 'N2OL'
517 !
518 ! :operatorSubType(2,i):Operator sub-type name.
519 ! Index (2,i) for sub-type to apply for stnid in element (1,i)
520 ! See related "obsoper@operatorCategory" automatically assigned
521 ! based on obs BUFR element and obsinfo_chm content.
522 !
523 ! Operator Sub-type name Description
524 ! Category
525 ! ============= =====================================================================
526 ! 'Interp' 'default' Piecewise linear interpolation (default)
527 ! 'wgtAvg' Piecewise weighted averaging interpolator
528 ! 'Surface' 'default' No special treatment (default)
529 ! 'Integ' 'default' Simple/basic vertical integration (default)
530 ! 'genOper' Weighted vertical integration - see 'genOper*' parameters
531 ! 'LayerAvg' 'default' Simple layer averaging (default)
532 ! 'genOper' Weighted vertical layer averaging - see 'genOper*' parameters
533 ! ====================================================================================
534 !
535 ! Notes:
536 !
537 ! - 'genOper' requires NAMBCHM namelist parameter settings
538 ! getPhysSpaceStats=.true. and
539 ! getPhysSpaceHCorrel=.true.
540 ! - Application of averaging kernels is directed only
541 ! by the content of the obsinfo_chm file 'SECTION III'
542 !
543 ! :storeOperators: Logical indicating if linear operators are stored for re-use in TL and AD calc.
544 ! If so, the linear operators will not be re-calculated at different iterations.
545 ! Not used when tropo_mode>=1
546 !
547 implicit none
548
549 ! Locals:
550 integer :: fnom, fclos
551 integer :: ierr, ios, nulnam, i
552 character(len=10) :: namfile
553
554 ! Namelist variables (local):
555 integer :: tropo_mode(0:oopc_constituentsSize) ! Special treatment for troposphere of total column obs
556 integer :: tropo_bound(0:oopc_constituentsSize) ! Indicate which column top value used for special treatment
557 real(8) :: tropo_column_top(0:oopc_constituentsSize) ! Default for column boundary (in Pa) of total column obs
558 logical :: storeOperators ! Choose to store linear operators for re-use in TL/AD
559 character(len=5) :: genOperConstraintType(0:oopc_constituentsSize) ! Strong constraint for generalized obs operator (see oopc_genOper)
560 real(8) :: genOperHCorrlenExpnt(0:oopc_constituentsSize) ! Exponent for horiz. correl. length weighting in oopc_genOper
561 real(8) :: genOperOmAStatsFactor(0:oopc_constituentsSize) ! Additional OmAStats normalization factor for oopc_genOper
562 integer :: obsdata_maxsize ! Max number of obs associated with ordered obs indices
563
564 external fnom,fclos
565
566 namelist /namchem/ assim_fam,assim_all,assim_num,assim_stnid,assim_varno, &
567 assim_nlev, assim_exclude_nflag,assim_exclude_flag, &
568 tropo_mode,tropo_bound,tropo_column_top,obsdata_maxsize, &
569 modelName,operatorSubType,storeOperators, &
570 genOperConstraintType,genOperHCorrlenExpnt,genOperOmAStatsFactor
571
572 ! Default NAMCHEM values
573
574 genOperConstraintType(:)='Trial'
575 genOperHCorrlenExpnt(:)=1.7
576 genOperOmAStatsFactor(:)=1.0
577
578 obsdata_maxsize=90000
579
580 storeOperators=.false.
581 assim_fam(:)=''
582 assim_fam(1)='CH'
583 assim_famNum=1
584 assim_all(:)=.true.
585 assim_num(:)=0
586 assim_stnid(:,:)='*********'
587 assim_varno(:,:)=0
588 assim_nlev(:,:)=0
589
590 assim_exclude_nflag(:)=1
591 assim_exclude_flag(:,1)=6 ! this is for the 'in reserve' bit for a BURP marker
592 assim_exclude_flag(:,2:)=0
593
594 tropo_mode(:) = 0
595 tropo_bound(:) = 0
596 tropo_column_top(:) = 0.0
597
598 operatorSubType(1,:) = ''
599 operatorSubType(2,:) = 'default'
600
601 ! Read from namelist file NAMCHEM
602
603 namfile=trim("flnml")
604 nulnam=0
605 ierr=fnom(nulnam,namfile,'R/O',0)
606
607 read(nulnam,nml=namchem,iostat=ios)
608 if (ios < -4 .or. ios > 0) then
609 call utl_abort('oopc_readNamchem: Error in reading NAMCHEM namelist. iostat = ' // trim(utl_str(ios)) )
610 else if (mmpi_myid == 0) then
611 write(*,nml=namchem)
612 end if
613
614 ierr=fclos(nulnam)
615
616 do i=size(assim_fam),1,-1
617 if (assim_fam(i) /= '') exit
618 end do
619 assim_famnum=i
620
621 oopc_genOperConstraintType(:) = genOperConstraintType(:)
622 oopc_genOperHCorrlenExpnt(:) = genOperHCorrlenExpnt(:)
623 oopc_genOperOmAStatsFactor(:) = genOperOmAStatsFactor(:)
624 oopc_tropo_mode(:) = tropo_mode(:)
625 oopc_tropo_bound(:) = tropo_bound(:)
626 oopc_tropo_column_top(:) = tropo_column_top(:)
627 oopc_storeOperators = storeOperators
628 oopc_obsdata_maxsize=obsdata_maxsize
629
630 end subroutine oopc_readNamchem
631
632 !--------------------------------------------------------------------------
633 !------- Routines related to reference or layer top & bottom levels -------
634
635 !--------------------------------------------------------------------------
636 ! oopc_readLevels
637 !--------------------------------------------------------------------------
638 subroutine oopc_readLevels
639 !
640 !:Purpose: To read and to store reference levels (where needed) or top- and
641 ! bottom-layer boundaries for CH sub-families.
642 !
643 implicit none
644
645 ! Locals:
646 integer :: fnom, fclos
647 integer :: ierr, jlev, jelm, nulstat, ios, isize, icount
648 logical :: LnewExists,newread
649 character (len=128) :: ligne
650 external :: fnom,fclos
651
652 ! Initialization
653
654 oopc_levels%n_stnid=0
655
656 inquire(file=trim(oopc_aux_filename),EXIST=LnewExists)
657 if (.not.LnewExists )then
658 write(*,*) '---------------------------------------------------------------'
659 write(*,*) 'oopc_readLevels: COULD NOT FIND AUXILIARY file ' // trim(oopc_aux_filename)
660 write(*,*) '---------------------------------------------------------------'
661 return
662 end if
663
664 ! Check for available layer info.
665
666 nulstat=0
667 ierr=fnom(nulstat,trim(oopc_aux_filename),'SEQ',0)
668 if ( ierr == 0 ) then
669 open(unit=nulstat, file=trim(oopc_aux_filename), status='OLD')
670 else
671 call utl_abort('oopc_readLevels: COULD NOT OPEN AUXILIARY file ' // trim(oopc_aux_filename))
672 end if
673
674 ios=0
675 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
676 do while (trim(adjustl(ligne(1:13))) /= 'SECTION II:')
677 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
678 end do
679
680 ! Read number of observation set sub-families (STNIDs and ...) and allocate space
681
682 read(nulstat,*,iostat=ios,err=10,end=10) oopc_levels%n_stnid
683 read(nulstat,*,iostat=ios,err=10,end=10) isize
684
685 allocate(oopc_levels%stnids(oopc_levels%n_stnid))
686 allocate(oopc_levels%vco(oopc_levels%n_stnid))
687 allocate(oopc_levels%source(oopc_levels%n_stnid),oopc_levels%ibegin(oopc_levels%n_stnid))
688 allocate(oopc_levels%element(oopc_levels%n_stnid),oopc_levels%n_lvl(oopc_levels%n_stnid))
689 allocate(oopc_levels%vlayertop(isize),oopc_levels%vlayerbottom(isize))
690 allocate(oopc_levels%n_col(oopc_levels%n_stnid))
691
692 oopc_levels%element(:)=0
693 oopc_levels%vco(:)=0
694 oopc_levels%source(:)=0
695 oopc_levels%n_lvl(:)=1
696 oopc_levels%n_col(:)=2
697
698 ! Begin reading for each sub-family
699 ! Important: Combination of STNID, BUFR element and number of vertical levels
700 ! to determine association to the observations.
701
702 icount=0
703 do jelm=1,oopc_levels%n_stnid
704 oopc_levels%ibegin(jelm)=icount+1
705
706 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
707
708 ! Read STNID (* is a wildcard)
709
710 read(nulstat,'(2X,A9)',iostat=ios,err=10,end=10) oopc_levels%stnids(jelm)
711
712 ! Read (1) Obs BUFR element.
713 ! (2) Vertical coord type (1, 2, or 3)
714 ! (3) Flag indication if EOR provided from this auxiliary file or
715 ! to be read from the observation file,
716 ! (4) Number of vertical levels (rows)
717 ! (5) Number of vertical level profiles/columns (optional; default is 2)
718 !
719 ! Important: Combination of STNID, BUFR element and number of vertical levels
720 ! to determine association to the observations.
721
722 newread=.true.
723 read(nulstat,*,iostat=ios,err=20,end=20) oopc_levels%element(jelm),oopc_levels%vco(jelm), &
724 oopc_levels%source(jelm),oopc_levels%n_lvl(jelm),oopc_levels%n_col(jelm)
725 newread=.false.
726 20 if (newread) then
727 backspace(nulstat,iostat=ios,err=10)
728 backspace(nulstat,iostat=ios,err=10)
729 read(nulstat,*,iostat=ios,err=10,end=10) oopc_levels%element(jelm),oopc_levels%vco(jelm), &
730 oopc_levels%source(jelm),oopc_levels%n_lvl(jelm)
731 end if
732
733 if (icount+oopc_levels%n_lvl(jelm) > isize) then
734 call utl_abort('oopc_readLevels: READING PROBLEM. ' // &
735 'Max array size exceeded: ' // &
736 trim(utl_str(isize)))
737 end if
738
739 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
740
741 if ( oopc_levels%n_lvl(jelm) > 0 ) then
742 if (oopc_levels%n_col(jelm) > 1) then
743
744 do jlev=1,oopc_levels%n_lvl(jelm)
745 icount=icount+1
746
747 ! Read top and bottom levels
748
749 read(nulstat,*,iostat=ios,err=10,end=10) &
750 oopc_levels%vlayertop(icount),oopc_levels%vlayerbottom(icount)
751 end do
752
753 else
754
755 do jlev=1,oopc_levels%n_lvl(jelm)
756 icount=icount+1
757
758 ! Read single level
759
760 read(nulstat,*,iostat=ios,err=10,end=10) &
761 oopc_levels%vlayertop(icount)
762
763 end do
764 oopc_levels%vlayerbottom(:)= oopc_levels%vlayertop(:)
765
766 end if
767 end if
768
769 end do
770
771 10 if (ios > 0) then
772 call utl_abort('oopc_readLevels: READING PROBLEM. ' // &
773 'File read error message number: ' // trim(utl_str(ios)))
774 end if
775
776 11 CLOSE(UNIT=nulstat)
777 ierr=fclos(nulstat)
778
779 end subroutine oopc_readLevels
780
781 !--------------------------------------------------------------------------
782 ! oopc_getLevels
783 !--------------------------------------------------------------------------
784 subroutine oopc_getLevels
785 !
786 !:Purpose: To return reference model levels or layer boundaries for an observation.
787 ! Combination of STNID, element variable number and number of vertical
788 ! levels to determine association to the observations. Default values
789 ! for top and bottom layers for total column measurements are to be
790 ! provided.
791 !
792 implicit none
793
794 ! Locals:
795 integer :: nlev ! number of levels in the observation
796 integer :: istnid,stnidIndex,startIndex,levelIndex,level
797 logical :: iset
798 real(8), allocatable :: levelsTop(:),levelsBot(:)
799
800 nlev=obsoper%nobslev
801
802 ! Find stnid with same number of vertical levels, and same BUFR element
803
804 istnid=0
805 obsoper%layerIdentified=.false.
806
807 do stnidIndex=1,oopc_levels%n_stnid
808
809 ! First compare STNID values allowing for * and blanks in
810 ! oopc_levels%stnids(JN) as wildcards
811 iset = utl_stnid_equal(oopc_levels%stnids(stnidIndex),obsoper%stnid)
812
813 ! Check if number of levels, code, and vertical coordinate type are equal.
814 ! If number of levels is one and no vertical coordinate provided for total column measurement (i.e. obsoper%vco == 4),
815 ! then check of vertical coordinate type is disregarded afterwards.
816 if (iset) then
817 if ( obsoper%varno == oopc_levels%element(stnidIndex) .and. &
818 (nlev == oopc_levels%n_lvl(stnidIndex) .or. &
819 (nlev == 1 .and. oopc_levels%n_lvl(stnidIndex) > 1) ) .and. &
820 (obsoper%vco == oopc_levels%vco(stnidIndex) .or. &
821 obsoper%vco == 4) ) then
822
823 istnid=stnidIndex
824 exit
825
826 end if
827 end if
828
829 end do
830
831 if (istnid == 0) then
832 ! If integrated layer information not found, if a total column measurement
833 ! set to defaults, else do nothing
834
835 if (bufr_IsIntegral(obsoper%varno) .and. nlev == 1) then
836 obsoper%layerIdentified=.true.
837 obsoper%vlayertop(1) = obsoper%pp(1)
838 obsoper%vlayerbottom(1) = obsoper%pp(obsoper%nmodlev)
839 end if
840
841 else
842 ! level or layer information has been found in auxiliary file
843
844 if (nlev == oopc_levels%n_lvl(stnidIndex) .and. nlev > 1) then
845 ! layer information has been found in auxiliary file
846 obsoper%layerIdentified=.true.
847 startIndex = oopc_levels%ibegin(stnidIndex)
848 obsoper%vlayertop(:) = oopc_levels%vlayertop(startIndex:startIndex+nlev-1)
849 obsoper%vlayerbottom(:) = oopc_levels%vlayerbottom(startIndex:startIndex+nlev-1)
850 else
851
852 ! Alternative level or layer information has been found
853 ! Must be pressure or sigma
854
855 if (oopc_levels%vco(stnidIndex) /= 2 &
856 .and. oopc_levels%vco(stnidIndex) /= 6) then
857
858 call utl_abort('oopc_getLevels: Cannot handle this vertical coordinate type')
859 end if
860
861 obsoper%nobslev=oopc_levels%n_lvl(stnidIndex)
862
863 if (oopc_levels%n_col(stnidIndex) > 1) obsoper%layerIdentified=.true.
864
865 startIndex = oopc_levels%ibegin(stnidIndex)
866
867 allocate(levelsTop(obsoper%nobslev))
868 allocate(levelsBot(obsoper%nobslev))
869 levelsTop(:) = oopc_levels%vlayertop(startIndex:startIndex+obsoper%nobslev-1)
870 levelsBot(:) = oopc_levels%vlayerbottom(startIndex:startIndex+obsoper%nobslev-1)
871 if (oopc_levels%vco(stnidIndex) == 6) then
872 levelsTop(:) = levelsTop(:)*obsoper%pp(obsoper%nmodlev)
873 levelsBot(:) = levelsBot(:)*obsoper%pp(obsoper%nmodlev)
874 end if
875
876 ! Ensure levels do not go below the surface
877
878 do levelIndex=1,obsoper%nobslev
879 if (levelsBot(levelIndex) >= obsoper%pp(obsoper%nmodlev)) then
880 obsoper%nobslev = levelIndex
881 levelsBot(levelIndex) = obsoper%pp(obsoper%nmodlev)
882 if (levelsTop(levelIndex) >= obsoper%pp(obsoper%nmodlev)) then
883 levelsTop(levelIndex) = obsoper%pp(obsoper%nmodlev)
884 end if
885 exit
886 end if
887 end do
888
889 ! Reset work array sizes and settings
890 deallocate(obsoper%vlayertop,obsoper%vlayerbottom,obsoper%obslev)
891 deallocate(obsoper%modlevindexTop,obsoper%modlevindexBot,obsoper%zh,obsoper%zhp)
892
893 allocate(obsoper%obslev(obsoper%nobslev)) ! Reference vertical levels
894 allocate(obsoper%vlayertop(obsoper%nobslev)) ! Layer tops for layer measurements
895 allocate(obsoper%vlayerbottom(obsoper%nobslev)) ! Layer bottoms for layer measurements
896 allocate(obsoper%modlevindexTop(obsoper%nobslev)) ! Index of highest model level (lowest index) involved with obs element
897 allocate(obsoper%modlevindexBot(obsoper%nobslev)) ! Index of lowest model level (highest index) involved with obs element
898 allocate(obsoper%zh(obsoper%nobslev,obsoper%nmodlev)) ! Local model operator H (excluding conversion constants and horizontal interpolation)
899 allocate(obsoper%zhp(obsoper%nobslev,obsoper%nmodlev)) ! Part of zh that excludes aspects related to vertical resolition
900
901 obsoper%vlayertop(:) = levelsTop(1:obsoper%nobslev)
902 obsoper%vlayerbottom(:) = levelsBot(1:obsoper%nobslev)
903 obsoper%obslev(:) = levelsTop(1:obsoper%nobslev)
904 obsoper%zh(:,:)=0.0d0
905 obsoper%zhp(:,:)=0.0d0
906
907 deallocate(levelsTop,levelsBot)
908
909 end if
910
911 obsoper%modlevindexTop(:)=1
912 startIndex=2
913 do level=1,obsoper%nobslev
914 if ( obsoper%vlayertop(level) < obsoper%pp(2)) then
915 obsoper%modlevindexTop(level)=1
916 else
917 do levelIndex=startIndex,obsoper%nmodlev
918 if (obsoper%vlayertop(level) == obsoper%pp(levelIndex)) then
919 obsoper%modlevindexTop(level)=levelIndex
920 exit
921 else if (obsoper%vlayertop(level) < obsoper%pp(levelIndex)) then
922 obsoper%modlevindexTop(level)=levelIndex-1
923 exit
924 end if
925 end do
926 if (level > 1) then
927 if (obsoper%modlevindexTop(level) < obsoper%modlevindexTop(level-1)) &
928 obsoper%modlevindexTop(level)=obsoper%modlevindexTop(level-1)
929 end if
930 startIndex=levelIndex
931 end if
932 end do
933
934
935 obsoper%modlevindexBot(:)=obsoper%modlevindexTop(:)+1
936 startIndex=obsoper%nmodlev-1
937 do level=obsoper%nobslev,1,-1
938 if ( obsoper%vlayertop(level) > obsoper%pp(obsoper%nmodlev-1)) then
939 obsoper%modlevindexTop(level)=obsoper%nmodlev
940 else
941 do levelIndex=startIndex,1,-1
942 if (obsoper%vlayerbottom(level) == obsoper%pp(levelIndex)) then
943 obsoper%modlevindexBot(level)=levelIndex
944 exit
945 else if (obsoper%vlayerbottom(level) > obsoper%pp(levelIndex)) then
946 obsoper%modlevindexBot(level)=levelIndex+1
947 exit
948 end if
949 end do
950 if (level < obsoper%nobslev) then
951 if (obsoper%modlevindexBot(level) > obsoper%modlevindexBot(level+1)) &
952 obsoper%modlevindexBot(level)=obsoper%modlevindexBot(level+1)
953 end if
954 startIndex=levelIndex
955 end if
956 end do
957
958 end if
959
960 end subroutine oopc_getLevels
961
962 !--------------------------------------------------------------------------
963 ! oopc_deallocLevels
964 !--------------------------------------------------------------------------
965 subroutine oopc_deallocLevels
966 !
967 !:Purpose: To deallocate temporary storage space used for layer info
968 !
969 implicit none
970
971 if (oopc_levels%n_stnid == 0) return
972
973 call oopc_deallocInfo(oopc_levels)
974
975 end subroutine oopc_deallocLevels
976
977 !--------------------------------------------------------------------------
978 !-------------- Routines related to averaging kernel matrices -------------
979
980 !--------------------------------------------------------------------------
981 ! oopc_readAvgkern
982 !--------------------------------------------------------------------------
983 subroutine oopc_readAvgkern
984 !
985 !:Purpose: To read averaging kernels from auxiliary file or observation file
986 !
987 implicit none
988
989 ! Locals:
990 integer, parameter :: ndim=2
991 integer :: istnid
992
993 ! read the averaging kernel information from the auxiliary file
994 call oopc_readAvgkernAuxfile
995
996 ! set size of observation file array
997 allocate(oopc_avgkern%obsSubSpace(oopc_avgkern%n_stnid))
998
999 ! read from observation file
1000 do istnid=1,oopc_avgkern%n_stnid
1001 if (oopc_avgkern%source(istnid) == 1) then
1002
1003 ! retrieve data from stats blocks (with bkstp=14 and block_type='DATA')
1004 oopc_avgkern%obsSubSpace(istnid) = obsf_obsSub_read('CH', &
1005 oopc_avgkern%stnids(istnid),bufr_avgkern,oopc_avgkern%n_lvl(istnid), &
1006 ndim, numColumns_opt=oopc_avgkern%n_col(istnid), &
1007 bkstp_opt=14, block_opt='DATA', match_nlev_opt=.true.)
1008
1009 end if
1010 end do
1011
1012 end subroutine oopc_readAvgkern
1013
1014 !--------------------------------------------------------------------------
1015 ! oopc_readAvgkernAuxfile
1016 !--------------------------------------------------------------------------
1017 subroutine oopc_readAvgkernAuxfile
1018 !
1019 !:Purpose: To read and to store averaging kernel matricesfor CH sub-families
1020 !
1021 !:Comments:
1022 ! - Currently implemented for only one latitude band
1023 !
1024
1025 implicit none
1026
1027 ! Locals:
1028 integer :: fnom, fclos
1029 integer :: ierr, jlev, jelm, nulstat, ios, isize, icount, iend
1030 logical :: LnewExists,newread
1031 character (len=128) :: ligne
1032 external :: fnom,fclos
1033
1034 ! Initialization
1035
1036 oopc_avgkern%n_stnid=0
1037
1038 inquire(file=trim(oopc_aux_filename),EXIST=LnewExists)
1039 if (.not.LnewExists )then
1040 write(*,*) '-------------------------------------------------------'
1041 write(*,*) 'oopc_readAvgkernAuxfile: COULD NOT FIND AUXILIARY file ' // trim(oopc_aux_filename)
1042 write(*,*) '-------------------------------------------------------'
1043 return
1044 end if
1045
1046 ! Check for available layer info.
1047
1048 nulstat=0
1049 ierr=fnom(nulstat,trim(oopc_aux_filename),'SEQ',0)
1050 if ( ierr == 0 ) then
1051 open(unit=nulstat, file=trim(oopc_aux_filename), status='OLD')
1052 else
1053 call utl_abort('oopc_readAvgkernAuxfile: COULD NOT OPEN AUXILIARY file ' // trim(oopc_aux_filename))
1054 end if
1055
1056 ios=0
1057 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
1058 do while (trim(adjustl(ligne(1:14))) /= 'SECTION III:')
1059 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
1060 end do
1061
1062 ! Read number of observation set sub-families (STNIDs and ...) and allocate space
1063
1064 read(nulstat,*,iostat=ios,err=10,end=10) oopc_avgkern%n_stnid
1065 read(nulstat,*,iostat=ios,err=10,end=10) isize
1066
1067 allocate(oopc_avgkern%stnids(oopc_avgkern%n_stnid))
1068 allocate(oopc_avgkern%source(oopc_avgkern%n_stnid))
1069 allocate(oopc_avgkern%ibegin(oopc_avgkern%n_stnid))
1070 allocate(oopc_avgkern%element(oopc_avgkern%n_stnid))
1071 allocate(oopc_avgkern%n_lvl(oopc_avgkern%n_stnid))
1072 allocate(oopc_avgkern%n_col(oopc_avgkern%n_stnid))
1073 allocate(oopc_avgkern%type(oopc_avgkern%n_stnid))
1074 allocate(oopc_avgkern%profElement(oopc_avgkern%n_stnid))
1075 allocate(oopc_avgkern%rak(isize))
1076
1077 oopc_avgkern%element(:)=0
1078 oopc_avgkern%source(:)=0
1079 oopc_avgkern%n_lvl(:)=1
1080 oopc_avgkern%n_col(:)=1
1081 oopc_avgkern%type(:) = 'default'
1082
1083 ! Begin reading for each sub-family
1084 ! Important: Combination of STNID, BUFR element and number of vertical levels
1085 ! to determine association to the observations.
1086
1087 icount=1
1088 STNIDLOOP: do jelm=1,oopc_avgkern%n_stnid
1089 oopc_avgkern%ibegin(jelm)=icount
1090
1091 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
1092
1093 ! Read STNID (* is a wildcard)
1094
1095 read(nulstat,'(2X,A9,1X,A15)',iostat=ios,err=10,end=10) oopc_avgkern%stnids(jelm),oopc_avgkern%type(jelm)
1096 if (trim(oopc_avgkern%type(jelm)) /= 'log') then
1097 oopc_avgkern%type(jelm) = 'default'
1098 end if
1099
1100 ! Read (1) Obs BUFR element.
1101 ! (2) Flag indication if avgkern provided from this auxiliary file or
1102 ! to be read from an observation file,
1103 ! (3) Number of obs in profile (number of rows)
1104 ! (4) Number of columns (optional read; number of obs operator vertical levels with or without the
1105 ! resultant a priori contribution (I-A)xa).
1106 ! If number of rows (n_lvl) equals 1, the actual levels for the columns need to be specified from
1107 ! section II (see routine oopc_readLevels)
1108 !
1109 ! Important: Combination of STNID, BUFR element and number of vertical levels
1110 ! to determine association to the observations.
1111
1112 newread=.true.
1113 read(nulstat,*,iostat=ios,err=20,end=20) oopc_avgkern%element(jelm), &
1114 oopc_avgkern%source(jelm),oopc_avgkern%n_lvl(jelm),oopc_avgkern%n_col(jelm),oopc_avgkern%profElement(jelm)
1115 newread=.false.
1116 20 if (newread) then
1117 backspace(nulstat,iostat=ios,err=10)
1118 backspace(nulstat,iostat=ios,err=10)
1119 read(nulstat,*,iostat=ios,err=10,end=10) oopc_avgkern%element(jelm), &
1120 oopc_avgkern%source(jelm),oopc_avgkern%n_lvl(jelm)
1121 oopc_avgkern%n_col(jelm)=oopc_avgkern%n_lvl(jelm)
1122 oopc_avgkern%profElement(jelm)= oopc_avgkern%element(jelm)
1123 end if
1124
1125 if (icount+oopc_avgkern%n_lvl(jelm) > isize) then
1126 call utl_abort('oopc_readAvgkernAuxfile: READING PROBLEM. ' // &
1127 'Max array size exceeded:' // trim(utl_str(isize)))
1128 end if
1129 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
1130
1131 ! disregard data section if values to be specified in BUFR file
1132 if (oopc_avgkern%source(jelm) == 1) cycle STNIDLOOP
1133
1134 if (oopc_avgkern%n_lvl(jelm) > 0 .and. oopc_avgkern%n_col(jelm) > 1) then
1135 do jlev=1,oopc_avgkern%n_lvl(jelm)
1136
1137 iend=icount+oopc_avgkern%n_col(jelm)-1
1138
1139 ! Read averaging kernel matrix
1140 read(nulstat,*,iostat=ios,err=10,end=10) oopc_avgkern%rak(icount:iend)
1141
1142 icount=iend+1
1143
1144 end do
1145 end if
1146
1147 end do STNIDLOOP
1148
1149 10 if (ios > 0) then
1150 call utl_abort('oopc_readAvgkernAuxfile: READING PROBLEM. ' // &
1151 'File read error message number: ' // trim(utl_str(ios)))
1152 end if
1153
1154 11 CLOSE(UNIT=nulstat)
1155 ierr=fclos(nulstat)
1156
1157 end subroutine oopc_readAvgkernAuxfile
1158
1159 !--------------------------------------------------------------------------
1160 ! oopc_deallocAvgkern
1161 !--------------------------------------------------------------------------
1162 subroutine oopc_deallocAvgkern
1163 !
1164 !:Purpose: To deallocate temporary storage space used for averaging kernels
1165 !
1166 implicit none
1167
1168 ! Locals:
1169 integer :: istnid
1170
1171 if (oopc_avgkern%n_stnid == 0) return
1172
1173 if (allocated(oopc_avgkern%obsSubSpace)) then
1174 do istnid=1,oopc_avgkern%n_stnid
1175 if (oopc_avgkern%source(istnid) == 1) then
1176 call oss_obsdata_dealloc(oopc_avgkern%obsSubSpace(istnid))
1177 end if
1178 end do
1179 deallocate(oopc_avgkern%obsSubSpace)
1180 end if
1181
1182 call oopc_deallocInfo(oopc_avgkern)
1183
1184 end subroutine oopc_deallocAvgkern
1185
1186 !--------------------------------------------------------------------------
1187 ! oopc_findAvgkern
1188 !--------------------------------------------------------------------------
1189 function oopc_findAvgkern(cstnid,varno,nlev) result(istnid)
1190 !
1191 !:Purpose: To find the averaging kernel for an observation if one is
1192 ! specified. Returns 0 if either not found or not specified.
1193 ! Combination of STNID, BUFR element and number of vertical levels
1194 ! to determine association to the observations.
1195 !
1196 implicit none
1197
1198 ! Arguments:
1199 character(len=12), intent(in) :: cstnid ! station id
1200 integer, intent(in) :: varno ! BUFR descriptor element
1201 integer, intent(in) :: nlev ! number of levels in the observation
1202 ! Result:
1203 integer :: istnid ! Index of averaging kernel in oopc_avgkern. Zero if not found.
1204
1205 ! Locals:
1206 integer :: stnidIndex
1207 logical :: iset
1208
1209 ! Find stnid with same number of vertical levels, and same BUFR element
1210
1211 istnid=0
1212
1213 do stnidIndex=1,oopc_avgkern%n_stnid
1214
1215 ! First compare STNID values allowing for * and blanks in
1216 ! oopc_avgkern%stnids(stnidIndex) as wildcards
1217 iset = utl_stnid_equal(oopc_avgkern%stnids(stnidIndex),CSTNID)
1218
1219 ! Check if number of levels and BUFR code are equal.
1220 if (iset) then
1221 if ( varno == oopc_avgkern%element(stnidIndex) ) then
1222 if (nlev < 1 .or. ( nlev == oopc_avgkern%n_lvl(stnidIndex) .or. &
1223 nlev == oopc_avgkern%n_col(stnidIndex) .or. &
1224 nlev == oopc_avgkern%n_col(stnidIndex)-1 .or. &
1225 nlev == oopc_avgkern%n_col(stnidIndex)-2) ) then
1226
1227 istnid=stnidIndex
1228 exit
1229 end if
1230 end if
1231 end if
1232
1233 end do
1234
1235 end function oopc_findAvgkern
1236
1237 !--------------------------------------------------------------------------
1238 ! oopc_getAvgkern
1239 !--------------------------------------------------------------------------
1240 subroutine oopc_getAvgkern(istnid,nlev,ncol,code,avg_kern)
1241 !
1242 !:Purpose: To return averaging kernel for an observation.
1243 !
1244 implicit none
1245
1246 ! Arguments:
1247 integer, intent(in) :: istnid ! index of averaging kernel in oopc_avgkern
1248 character(len=*), intent(in) :: code ! measurement identifier
1249 integer, intent(in) :: nlev ! number of observation levels
1250 integer, intent(in) :: ncol ! number of columns for avg kernel info (without the a priori contribution)
1251 real(8), intent(out) :: avg_kern(nlev,ncol+2) ! averaging kernel, plus possible a priori contribution and integration weights
1252
1253 ! Locals:
1254 integer :: startIndex,endIndex
1255
1256 if (istnid > 0 .and. istnid <= oopc_avgkern%n_stnid) then
1257
1258 if (oopc_avgkern%source(istnid) == 0) then
1259 ! Check number of columns
1260 if (ncol < oopc_avgkern%n_col(istnid) .or. &
1261 ncol+2 > oopc_avgkern%n_col(istnid) ) then
1262
1263 call utl_abort('oopc_getAvgkern: Inconsistency ' // &
1264 'in avg kern size for ' // oopc_avgkern%stnids(istnid) )
1265 end if
1266
1267 ! get averaging kernel from auxiliary file
1268 startIndex = oopc_avgkern%ibegin(istnid)
1269 endIndex = nlev*(startIndex+oopc_avgkern%n_col(istnid)-1)
1270 avg_kern = RESHAPE(oopc_avgkern%rak(startIndex:endIndex), &
1271 (/nlev,oopc_avgkern%n_col(istnid)/),ORDER =(/2,1/))
1272 else
1273 ! get averaging kernel from observation file
1274 avg_kern(1:nlev,1:oopc_avgkern%n_col(istnid)) = &
1275 oss_obsdata_get_array2d(oopc_avgkern%obsSubSpace(istnid),code)
1276 end if
1277
1278 else
1279 call utl_abort("oopc_getAvgkern: Invalid station ID index.")
1280 end if
1281
1282 end subroutine oopc_getAvgkern
1283
1284 !----------------------------------- Misc ---------------------------------
1285 !--------------------------------------------------------------------------
1286
1287 !--------------------------------------------------------------------------
1288 ! oopc_dealocInfo
1289 !--------------------------------------------------------------------------
1290 subroutine oopc_deallocInfo(info)
1291 !
1292 !:Purpose: To deallocate struct_oopc_info instance
1293 !
1294 implicit none
1295
1296 ! Arguments:
1297 type(struct_oopc_info), intent(inout) :: info
1298
1299 if (allocated(info%stnids)) deallocate(info%stnids)
1300 if (allocated(info%element)) deallocate(info%element)
1301 if (allocated(info%source)) deallocate(info%source)
1302 if (allocated(info%vco)) deallocate(info%vco)
1303 if (allocated(info%n_lat)) deallocate(info%n_lat)
1304 if (allocated(info%ibegin)) deallocate(info%ibegin)
1305 if (allocated(info%n_lvl)) deallocate(info%n_lvl)
1306 if (allocated(info%rak)) deallocate(info%rak)
1307 if (allocated(info%vlayertop)) deallocate(info%vlayertop)
1308 if (allocated(info%vlayerbottom)) deallocate(info%vlayerbottom)
1309 if (allocated(info%lat)) deallocate(info%lat)
1310
1311 end subroutine oopc_deallocInfo
1312
1313 !--------------------------------------------------------------------------
1314 ! oopc_diagnOnly
1315 !--------------------------------------------------------------------------
1316 logical function oopc_diagnOnly(cfamName,cstnid,varno,nobslev,flag)
1317 !
1318 !:Purpose: To identify whether or not the obs set identified by the
1319 ! combination of (cstnid,varno,nobslev) will be assimilated or
1320 ! else used for independent verifications after
1321 ! assimilation/minimization
1322 !
1323 implicit none
1324
1325 ! Arguments:
1326 character(len=*), intent(in) :: cfamName ! Family name
1327 character(len=*), intent(in) :: cstnid ! Input station id
1328 integer, intent(in) :: varno ! Obs BUFR number
1329 integer, intent(in) :: nobslev ! Number of levels
1330 integer, intent(in) :: flag ! observation integer flag
1331
1332 ! Locals:
1333 integer :: i,elemId,ifam
1334
1335 ifam=0
1336 if (assim_famNum > 0) then
1337 do i=1,assim_famNum
1338 if (assim_fam(i) == cfamName) then
1339 ifam=i
1340 exit
1341 end if
1342 end do
1343 end if
1344
1345 if (ifam == 0) then
1346 ! assimilate all observations
1347 oopc_diagnOnly = .false.
1348 return
1349 end if
1350
1351 if (assim_all(ifam)) then
1352 ! assimilate all observations
1353 oopc_diagnOnly = .false.
1354 else if (assim_num(ifam) <= 0) then
1355 ! assimilate no observations
1356 oopc_diagnOnly = .true.
1357 else if (assim_num(ifam) > 0) then
1358 ! check if this observation is listed in the assim_* arrays
1359 elemId=0
1360 do i=1,assim_num(ifam)
1361 if (utl_stnid_equal(trim(assim_stnid(ifam,i)),trim(cstnid))) then
1362 if (assim_varno(ifam,i) == 0 .or. assim_varno(ifam,i) == varno) then
1363 if (assim_nlev(ifam,i) == 0 .or. (nobslev == 1 .and. &
1364 assim_nlev(ifam,i) == 1) .or. &
1365 (nobslev > 1 .and. assim_nlev(ifam,i) > 1)) then
1366
1367 elemId=i
1368 exit
1369 end if
1370 end if
1371 end if
1372 end do
1373 oopc_diagnOnly = elemId == 0
1374 end if
1375
1376 if (oopc_diagnOnly) return
1377
1378 ! check if the observation integer flag has a bit marked by assim_exclude_flag (same flagging as in filt_suprep)
1379 if (assim_exclude_nflag(ifam) > 0) then
1380 do i=1,assim_exclude_nflag(ifam)
1381 if (btest(flag, 13 - assim_exclude_flag(ifam,i) )) then
1382 oopc_diagnOnly = .true.
1383 return
1384 end if
1385 end do
1386 end if
1387
1388 end function oopc_diagnOnly
1389
1390 !--------------------------------------------------------------------------
1391 ! oopc_checkType
1392 !--------------------------------------------------------------------------
1393 function oopc_checkType(StnidSet,TypeSet,stnid,type) result(sameType)
1394 !
1395 !:Purpose: To determine if specified combination of (stnid,type) found
1396 ! in (StnidSet,TypeSet).
1397 !
1398 implicit none
1399
1400 ! Arguments:
1401 character(len=*), intent(in) :: StnidSet(:)
1402 character(len=*), intent(in) :: TypeSet(:)
1403 character(len=*), intent(in) :: stnid
1404 character(len=*), intent(in) :: type
1405 ! Result:
1406 logical :: sametype
1407
1408 ! Locals:
1409 integer :: stnidIndex
1410
1411 sameType = .false.
1412
1413 do stnidIndex=1,size(StnidSet)
1414
1415 if ( trim(TypeSet(stnidIndex)) == '' .or. &
1416 trim(StnidSet(stnidIndex)) == '' ) exit
1417
1418 if (utl_stnid_equal(StnidSet(stnidIndex),stnid)) then
1419 if ( trim(StnidSet(stnidIndex)) == trim(stnid) ) then
1420 if ( trim(TypeSet(stnidIndex)) == trim(type) ) then
1421 sameType=.true.
1422 else
1423 sameType=.false.
1424 end if
1425 exit
1426 else
1427 if ( trim(TypeSet(stnidIndex)) == trim(type) ) sameType=.true.
1428 end if
1429 end if
1430
1431 end do
1432
1433 if ( .not.sameType .and. trim(type) == 'default' ) sameType = .true.
1434
1435 end function oopc_checkType
1436
1437 !--------------------------------------------------------------------------
1438 ! oopc_getType
1439 !--------------------------------------------------------------------------
1440 function oopc_getType(StnidSet,TypeSet,stnid) result(type)
1441 !
1442 !:Purpose: To determine "type" for specified "stnid" found in "(StnidSet,Typesef)".
1443 !
1444 implicit none
1445
1446 ! Arguments:
1447 character(len=*), intent(in) :: StnidSet(:)
1448 character(len=*), intent(in) :: TypeSet(:)
1449 character(len=*), intent(in) :: stnid
1450 ! Result:
1451 character(len=15) :: type
1452
1453 ! Locals:
1454 integer :: stnidIndex
1455
1456 type = 'default'
1457
1458 do stnidIndex=1,size(StnidSet)
1459
1460 if ( trim(TypeSet(stnidIndex)) == '' .or. &
1461 trim(StnidSet(stnidIndex)) == '' ) exit
1462
1463 if (utl_stnid_equal(StnidSet(stnidIndex),stnid)) then
1464 type=trim(TypeSet(stnidIndex))
1465 if ( trim(StnidSet(stnidIndex)) == trim(stnid) ) exit
1466 end if
1467
1468 end do
1469
1470 end function oopc_getType
1471
1472 !--------------------------------------------------------------------------
1473 !------------------ Routines associated to oopc_efftemp --------------------
1474
1475 !--------------------------------------------------------------------------
1476 ! oopc_addEfftempObsdata
1477 !--------------------------------------------------------------------------
1478 subroutine oopc_addEfftempObsdata(code,temp_eff)
1479 !
1480 !:Purpose: To add effective temperature value to its obsdata object
1481 !
1482 implicit none
1483
1484 ! Arguments:
1485 character(len=*), intent(in) :: code ! unique identifying code
1486 real(8), intent(in) :: temp_eff(:) ! effective temperature
1487
1488 call oss_obsdata_add_data1d(oopc_efftemp,temp_eff,code,oopc_obsdata_maxsize)
1489
1490 end subroutine oopc_addEfftempObsdata
1491
1492 !--------------------------------------------------------------------------
1493 ! oopc_addEfftempObsfile
1494 !--------------------------------------------------------------------------
1495 subroutine oopc_addEfftempObsfile()
1496 !
1497 !:Purpose: To add effective temperatures in obs file.
1498 !
1499 implicit none
1500
1501 ! Locals:
1502 integer :: nrep_modified,varno(1)
1503
1504 ! If needed, add effective temperature values in obs file for total column measurements
1505
1506 call oss_obsdata_MPIallgather(oopc_efftemp)
1507
1508 if (oopc_efftemp%nrep > 0) then
1509 varno(1)=12001
1510 nrep_modified = obsf_obsSub_update(oopc_efftemp,'CH', &
1511 varno(1:max(1,oopc_efftemp%dim2)),bkstp_opt=0, &
1512 block_opt='INFO',multi_opt='UNI')
1513 write(*,*) 'oopc_addEfftempObsfile: Added ',nrep_modified, &
1514 ' effective temperature values in the obs file.'
1515 end if
1516
1517 end subroutine oopc_addEfftempObsfile
1518
1519 !-------------- CONTROL ROUTINES FOR OBSERVATION OPERATORS ----------------
1520 !--------------------------------------------------------------------------
1521
1522 !--------------------------------------------------------------------------
1523 ! oopc_CHobsoperators
1524 !--------------------------------------------------------------------------
1525 subroutine oopc_CHobsoperators(columnTrl, obsSpaceData, mode, &
1526 columnAnlInc_opt, jobs_opt, destObsColumn_opt)
1527 !
1528 !:Purpose: To apply the observation operators for chemical constituents.
1529 ! Mode of operator set by mode (see also kmode below).
1530 !
1531 !:Comments:
1532 ! - See type struct_oopc_obsoperators for description of obsoper elements.
1533 ! - Currently can only handle the case when nlev_bkgrnd == nlev_inc
1534 !
1535 !:Arguments:
1536 ! :columnTrl: Column of x_background interpolated to observation
1537 ! location. Can have the same vertical levels as the
1538 ! trial field (columnhr) or as the increment field
1539 ! (columng)
1540 ! :mode: (kmode retained internally following the switch to mode as input)
1541 ! +------+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
1542 ! | mode |kmode| Mode of | Results |
1543 ! | | | Observation Operator | |
1544 ! +=======+=====+=======================+===================================+
1545 ! | nl | 0 |for general simulation |OmP and total Jo(x_background) |
1546 ! | | |operator |for CH. OmP saved in OBS_OMP of |
1547 ! | | |(non-linear and linear)|obsSpaceData |
1548 ! +-------+-----+-----------------------+-----------------------------------+
1549 ! | HBHT | 1 |for identification or |background error standard dev. |
1550 ! | | |determination of |in observation space saved in |
1551 ! | | |of sqrt(diag(H*B*H^T). |in OBS_HPHT of obsSpaceData |
1552 ! | | |Depends on the presence|if OmP error std dev not initially
1553 ! | | |of OmP error std dev |available in OBS_OMPE |
1554 ! +-------+-----+-----------------------+-----------------------------------+
1555 ! | tl | 2 |for tangent linear |Hdx saved in OBS_WORK of |
1556 ! | | |operator |obsSpaceData |
1557 ! +------+-----+-----------------------+-----------------------------------+
1558 ! |adjoint| 3 |for adjoint of tangent |H^T * R^-1 (OmP-Hdx) in |
1559 ! | | |linear operator |columnAnlInc_opt |
1560 ! +-------+-----+-----------------------+-----------------------------------+
1561 !
1562 ! :columnAnlInc_opt: Optional argument for input/output of column of
1563 ! increment (column). For kmode=2, used as input for
1564 ! increment H_horiz dx interpolated to observation
1565 ! location. For kmode=3, used as output for H^T * R^-1
1566 ! (OmP-Hdx). Required for kmode=2,3.
1567 !
1568 !
1569 ! :jobs_opt: Optional output of total Jo(x_background) for chemical
1570 ! constituents. Required for kmode=0 and not provided
1571 ! otherwise.
1572 !
1573 !
1574
1575
1576 ! More Comments (not rendered in Sphinx):
1577 ! - Two equivalent methods for looping over a report body.
1578 !
1579 ! Method 1:
1580 !
1581 ! call obs_set_current_body_list(obsSpaceData,headerIndex)
1582 ! BODY: do
1583 !
1584 ! bodyIndex = obs_getBodyIndex(obsSpaceData)
1585 ! if (bodyIndex < 0) exit BODY1
1586 !
1587 ! ... obs_bodyElem_r(obsSpaceData, ... ,bodyIndex)
1588 !
1589 ! end do BODY
1590 !
1591 ! Method 2:
1592 !
1593 ! bodyIndex_start = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
1594 ! bodyIndex_end = bodyIndex_start+obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)-1
1595 ! do bodyIndex=bodyIndex_start,bodyIndex_end
1596 ! ... obs_bodyElem_r(obsSpaceData, ... ,bodyIndex)
1597 ! end do
1598 !
1599 implicit none
1600
1601 ! Arguments:
1602 type(struct_columnData), intent(in) :: columnTrl
1603 type(struct_obs), intent(inout) :: obsSpaceData ! Observation-space data object
1604 character(len=*), intent(in) :: mode
1605 type(struct_columnData), optional, intent(in) :: columnAnlInc_opt
1606 real(8), optional, intent(out) :: jobs_opt
1607 integer, optional, intent(in) :: destObsColumn_opt
1608
1609 ! Locals:
1610 real(8) :: zomp,zinc,zoer,zhbht
1611 integer :: kmode ! Mode of observation operator
1612 integer, external :: fclos
1613 integer :: headerIndex,bodyIndex,bodyIndex_start,bodyIndex_end
1614 integer :: icodtyp,obslevIndex,nobslev,varno,maxnumHeaders,headerCount
1615 integer :: destObsColumn
1616 character(len=12) :: stnid
1617 integer, allocatable :: iass(:),flag(:)
1618 logical, allocatable :: process_obs(:)
1619 real(8), allocatable :: obs_col(:)
1620 real(8), pointer :: col(:),model_col(:)
1621 integer :: nlev_bkgrnd,nlev_inc,modlevIndex
1622 character(len=2), parameter :: varLevel = 'TH'
1623
1624 if ( mode == 'nl' ) then
1625 kmode = 0
1626 else if ( mode == 'HBHT' ) then
1627 kmode = 1
1628 else if ( mode == 'tl' ) then
1629 kmode = 2
1630 else if ( mode == 'adjoint' ) then
1631 kmode = 3
1632 end if
1633
1634 ! Apply setup on first call
1635 if (.not.initializedChem) then
1636 call oopc_setupCH(kmode)
1637 initializedChem = .true.
1638 else if ( mode == 'HBHT' .and. .not.bgStats%initialized ) then
1639 call bcsc_getCovarCH(bgStats)
1640 end if
1641
1642 if ((kmode == 2 .or. kmode == 3) .and. (.not.present(columnAnlInc_opt))) then
1643 call utl_abort("oopc_CHobsoperators: columnAnlInc_opt must " // &
1644 "be specified for kmode = " // utl_str(kmode))
1645 end if
1646
1647 ! Initializations
1648
1649 if ( present(destObsColumn_opt) ) then
1650 destObsColumn = destObsColumn_opt
1651 else
1652 destObsColumn = obs_omp
1653 end if
1654
1655 if (present(jobs_opt)) jobs_opt = 0.d0
1656
1657 nlev_bkgrnd = col_getNumLev(columnTrl,varLevel)
1658
1659 ! Allocate memory for model_col. Not necessary for kmode=0 since model_col points to obsoper%trial.
1660 select case(kmode)
1661 case(2)
1662 nlev_inc = col_getNumLev(columnAnlInc_opt,varLevel)
1663 allocate(model_col(nlev_inc))
1664 if (nlev_inc /= nlev_bkgrnd) then
1665 write(*,*) "oopc_CHobsoperators: nlev_inc ", &
1666 "and nlev_bkgrnd not the same: ",nlev_inc, nlev_bkgrnd
1667 end if
1668 case(1,3)
1669 allocate(model_col(nlev_bkgrnd))
1670 end select
1671
1672 ! Allocations outside oopc_obsoperInit since this can be done outside the HEADER loop.
1673 ! See oopc_obsoperInit for assignment of array content.
1674
1675 ! Model obs background, height, TT, and HU profiles.
1676 allocate(obsoper%trial(nlev_bkgrnd),obsoper%height(nlev_bkgrnd),obsoper%tt(nlev_bkgrnd),obsoper%hu(nlev_bkgrnd))
1677 ! Model PP and pressure model layer boundaries taken as the middle between model levels.
1678 allocate(obsoper%pp(nlev_bkgrnd))
1679
1680 ! Determine number of obs
1681 maxnumHeaders=0
1682 call obs_set_current_header_list(obsSpaceData,'CH')
1683 do
1684 maxnumHeaders=maxnumHeaders+1
1685 if (obs_getHeaderIndex(obsSpaceData) < 0) exit
1686 end do
1687
1688 ! Loop over all header indices of the 'CH' family:
1689
1690 headerCount=0
1691 call obs_set_current_header_list(obsSpaceData,'CH')
1692 HEADER: do
1693
1694 headerCount = headerCount+1
1695 headerIndex = obs_getHeaderIndex(obsSpaceData)
1696 if (headerIndex < 0) exit HEADER
1697
1698 icodtyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
1699 if (icodtyp /= codtyp_get_codtyp('CHEMREMOTE') .and. &
1700 icodtyp /= codtyp_get_codtyp('CHEMINSITU')) cycle HEADER
1701
1702 stnid = obs_elem_c(obsSpaceData,'STID',headerIndex)
1703 bodyIndex_start = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
1704 bodyIndex_end = bodyIndex_start+obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)-1
1705
1706 ! Set number of obs profile elements by removing count of BUFR_SCALE_EXPONENT elements
1707 nobslev = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)
1708 do bodyIndex=bodyIndex_start,bodyIndex_end
1709 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == &
1710 BUFR_SCALE_EXPONENT) then
1711 nobslev = nobslev-1
1712 end if
1713 end do
1714
1715 ! varno is expected to be the same for all profile points where OBS_VNM value /= BUFR_SCALE_EXPONENT
1716 do bodyIndex=bodyIndex_start,bodyIndex_end
1717 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) &
1718 /= BUFR_SCALE_EXPONENT) then
1719
1720 varno = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
1721 exit
1722 end if
1723 end do
1724
1725 ! Allocate memory for remaining profile data not in obsoper
1726 allocate(obs_col(nobslev),iass(nobslev),process_obs(nobslev),flag(nobslev))
1727
1728 if (allocated(obsoper%success)) deallocate(obsoper%success,obsoper%ixtr)
1729 allocate(obsoper%success(nobslev),obsoper%ixtr(nobslev))
1730
1731 ! Check to see if background error variances available
1732 if (kmode == 1) then
1733 process_obs(:) = bcsc_StatsExistForVarName(vnl_varnameFromVarnum(varno, &
1734 obs_headElem_i(obsSpaceData,OBS_CHM,headerIndex), &
1735 modelName))
1736 end if
1737
1738 ! Prepare for checking if any processing is needed according to initial flag values
1739 obslevIndex=0
1740
1741 do bodyIndex=bodyIndex_start,bodyIndex_end
1742 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) &
1743 /= BUFR_SCALE_EXPONENT) then
1744
1745 obslevIndex=obslevIndex+1
1746
1747 obsoper%ixtr(obslevIndex) = &
1748 obs_bodyElem_i(obsSpaceData,OBS_XTR,bodyIndex) ! indicates if obs extends outside model profile vertical range
1749 iass(obslevIndex) = obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) ! indicates if obs is to be assimilated
1750 flag(obslevIndex) = obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex) ! observation integer flag
1751
1752 ! Indicates if this obs should be processed by oopc_obsoperators
1753 if (kmode == 1) then
1754 process_obs(obslevIndex) = obsoper%ixtr(obslevIndex) == 0 &
1755 .and. iass(obslevIndex) == obs_assimilated .and. &
1756 process_obs(obslevIndex)
1757 else
1758 process_obs(obslevIndex) = obsoper%ixtr(obslevIndex) == 0 &
1759 .and. iass(obslevIndex) == obs_assimilated
1760 end if
1761
1762 end if
1763 end do
1764
1765 ! Initialize processing success flag
1766
1767 obsoper%success(1:nobslev) = process_obs(1:nobslev)
1768
1769 if (all(.not.process_obs)) then
1770
1771 ! All observations in the profile flagged so can skip obs operator for current measurement
1772
1773 if (kmode == 3) then
1774 model_col(:) = 0.0D0
1775 obsoper%varName = vnl_varnameFromVarnum(varno, &
1776 obs_headElem_i(obsSpaceData,OBS_CHM,headerIndex),modelName)
1777 end if
1778
1779 else
1780
1781 if ( kmode == 1 ) then
1782 ! Check if "OmP error std dev" is already available
1783 call obs_set_current_body_list(obsSpaceData,headerIndex)
1784 obslevIndex=0
1785 BODYINDEX1: do bodyIndex=bodyIndex_start,bodyIndex_end
1786
1787 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == &
1788 BUFR_SCALE_EXPONENT) cycle
1789
1790 obslevIndex = obslevIndex + 1
1791 if ( process_obs(obslevIndex) ) then
1792 if (obs_bodyElem_r(obsSpaceData,OBS_OMPE,bodyIndex) > 0.0d0 ) then
1793 ! "OmP error std dev" is already available for this measurement.
1794 ! Go to the next measurement.
1795
1796 ! TEMPORARY: First, estimate OBS_HPHT for storage in output files (in the event it is needed externally for
1797 ! other purposes (e.g. total column ozone bias correction and corresponding re-doing for marker settings)
1798 if ( obs_bodyElem_r(obsSpaceData,OBS_OMPE,bodyIndex) > &
1799 1.1d0*obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex) ) then
1800
1801 zhbht = sqrt(obs_bodyElem_r(obsSpaceData,OBS_OMPE,bodyIndex)**2 &
1802 -obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)**2)
1803 call obs_bodySet_r(obsSpaceData,OBS_HPHT,bodyIndex,zhbht)
1804 else
1805 call obs_bodySet_r(obsSpaceData,OBS_HPHT,bodyIndex, &
1806 0.5d0*obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex))
1807 end if
1808
1809 deallocate(process_obs,obsoper%success,obsoper%ixtr,iass,obs_col,flag)
1810
1811 cycle HEADER
1812 else
1813 ! Proceed with the calc of sqrt(diag(HBHT))
1814 exit BODYINDEX1
1815 end if
1816 end if
1817 end do BODYINDEX1
1818 end if
1819
1820 ! Initialize obsoper variables and allocate arrays
1821 call oopc_obsoperInit(obsSpaceData,headerIndex,columnTrl,nlev_bkgrnd,nobslev,kmode,varno,stnid)
1822
1823 ! Initialize model_col, dependent on kmode. Used for input for kmode=0,2, output for kmode=3.
1824 ! model_col represents for kmode 0) the horizontally interpolated background H_horiz(x_b)
1825 ! 1) not used
1826 ! 2) the analysis increment H_horiz dx
1827 ! 3) the result of applying the adjoint of H_vert
1828
1829 select case(kmode)
1830 case(0)
1831 model_col => obsoper%trial
1832 case(2)
1833 do modlevIndex=1,nlev_inc
1834 model_col(modlevIndex) = col_getElem(columnAnlInc_opt,modlevIndex,headerIndex,obsoper%varName)
1835 end do
1836 case(1,3)
1837 model_col(:) = 0.0D0
1838 end select
1839
1840 ! Loop over all body indices (profile elements) to acquire remaining data
1841
1842 obslevIndex=0
1843 call obs_set_current_body_list(obsSpaceData,headerIndex)
1844 BODY1: do
1845
1846 bodyIndex = obs_getBodyIndex(obsSpaceData)
1847 if (bodyIndex < 0) exit BODY1
1848
1849 ! Get position in profile and skip over BUFR_SCALE_EXPONENT elements
1850
1851 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) /= &
1852 BUFR_SCALE_EXPONENT) then
1853
1854 obslevIndex=obslevIndex+1
1855 else
1856 cycle BODY1
1857 end if
1858
1859 ! Get vertical coordinate data. Valid for point data values in profiles.
1860 ! For layer data values, vertical coordinate data will instead be assigned within oopc_obsoperators.
1861
1862 obsoper%obslev(obslevIndex) = &
1863 obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1864
1865 ! Get normalized increment
1866 if (kmode == 3) then
1867 if (iass(obslevIndex) == 1) then
1868 obs_col(obslevIndex) = obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
1869 else
1870 obs_col(obslevIndex) = 0.0D0
1871 end if
1872 end if
1873
1874 ! Get obs estimate from trial field (if kmode>1)
1875
1876 if (kmode <= 1) then
1877 obsoper%obsSpaceTrial(obslevIndex) = 0.0
1878 else
1879 ! Store for use by TL and AD of non-linear operators
1880 obsoper%obsSpaceTrial(obslevIndex) = &
1881 obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex) &
1882 - obs_bodyElem_r(obsSpaceData,OBS_OMP,bodyIndex)
1883 end if
1884
1885 end do BODY1
1886
1887 ! Apply observation operator. model_col,obs_col are the inputs/outputs in model,observation
1888 ! space, respectively. Other required inputs are in obsoper. Input/output is as follows:
1889 !
1890 ! kmode model_col obs_col
1891 ! ----- --------- -------
1892 ! 0 in out
1893 ! 1 not used out
1894 ! 2 in out
1895 ! 3 out in
1896
1897 call oopc_obsoperators(model_col,obs_col,maxnumHeaders,headerCount,kmode)
1898
1899 end if
1900
1901 ! Output results
1902
1903 if (kmode == 3) then
1904 ! Store H^T * R^-1 (OmP-Hdx) in columnInc
1905
1906 col => col_getColumn(columnAnlInc_opt,headerIndex,obsoper%varName)
1907 col(1:nlev_bkgrnd) = model_col(1:nlev_bkgrnd)
1908
1909 else
1910 ! Store results in obsSpaceData
1911
1912 obslevIndex=0
1913 call obs_set_current_body_list(obsSpaceData,headerIndex)
1914 BODY2: do
1915 bodyIndex = obs_getBodyIndex(obsSpaceData)
1916 if (bodyIndex < 0) exit BODY2
1917
1918 ! Get position in profile and skip over BUFR_SCALE_EXPONENT elements
1919 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) /= &
1920 BUFR_SCALE_EXPONENT) then
1921
1922 obslevIndex=obslevIndex+1
1923 else
1924 cycle BODY2
1925 end if
1926
1927 ! Check for success in calculations
1928 if (process_obs(obslevIndex) .and. .not.obsoper%success(obslevIndex)) then
1929 ! Observation was flagged within this call of oopc_CHobsoperators
1930 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
1931 call obs_bodySet_r(obsSpaceData,OBS_OMP,bodyIndex,0.0D0)
1932 call obs_bodySet_r(obsSpaceData,OBS_OMA,bodyIndex,0.0D0)
1933 call obs_bodySet_r(obsSpaceData,OBS_HPHT,bodyIndex,0.0D0)
1934 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex,0.0D0)
1935 call obs_bodySet_i(obsSpaceData,OBS_FLG,bodyIndex, &
1936 ibset(obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex),9) )
1937 cycle BODY2
1938 else if (iass(obslevIndex) == 0) then
1939 ! Observation was flagged previous to this call of oopc_CHobsoperators
1940 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
1941 call obs_bodySet_r(obsSpaceData,OBS_OMP,bodyIndex,0.0D0)
1942 call obs_bodySet_r(obsSpaceData,OBS_OMA,bodyIndex,0.0D0)
1943 call obs_bodySet_r(obsSpaceData,OBS_HPHT,bodyIndex,0.0D0)
1944 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex,0.0D0)
1945 cycle BODY2
1946 else if (.not.process_obs(obslevIndex) .and. kmode == 1) then
1947 call obs_bodySet_r(obsSpaceData,OBS_HPHT,bodyIndex,0.0D0)
1948 cycle BODY2
1949 end if
1950
1951 ! Store result in appropriate location in obsSpaceData
1952 select case(kmode)
1953 case(0)
1954
1955 ! Store OmP in OBS_OMP and add to Jo(x_background) of CH.
1956
1957 zomp = obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex) - obs_col(obslevIndex)
1958 call obs_bodySet_r(obsSpaceData,destObsColumn,bodyIndex,zomp)
1959
1960 if (oopc_diagnOnly('CH',stnid,varno,nobslev,flag(obslevIndex))) then
1961 ! Observation is for diagnostics and is not to be assimilated
1962 call obs_bodySet_i(obsSpaceData,OBS_ASS,bodyIndex,obs_notAssimilated)
1963 else if (present(jobs_opt) .and. iass(obslevIndex) == 1) then
1964 ! Add to Jo contribution (factor of 0.5 to be applied outside report loop)
1965 zinc = zomp/obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)
1966 jobs_opt = jobs_opt + zinc**2
1967 end if
1968
1969 case(1)
1970
1971 ! Background error standard deviations in
1972 ! observation space, sqrt(diag(H*B_static*H^T)),
1973 ! saved in OBS_HPHT of obsSpaceData.
1974 ! Resulting OmP error std dev estimate saved in OBS_OMPE
1975
1976 call obs_bodySet_r(obsSpaceData,OBS_HPHT,bodyIndex,obs_col(obslevIndex))
1977
1978 zoer = obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)
1979 zomp = sqrt(obs_col(obslevIndex)*obs_col(obslevIndex) + zoer*zoer)
1980 call obs_bodySet_r(obsSpaceData,OBS_OMPE,bodyIndex,zomp)
1981
1982 case(2)
1983
1984 ! Store Hdx in OBS_WORK of obsSpaceData
1985 call obs_bodySet_r(obsSpaceData,OBS_WORK,bodyIndex,obs_col(obslevIndex))
1986 end select
1987
1988 end do BODY2
1989
1990 end if
1991
1992 ! Deallocate profile data
1993 deallocate(process_obs,iass,obs_col,flag)
1994 call oopc_obsoper_dealloc
1995
1996 end do HEADER
1997
1998 deallocate(obsoper%trial,obsoper%pp,obsoper%tt)
1999 deallocate(obsoper%height,obsoper%hu)
2000 if (kmode.ne.0) deallocate(model_col)
2001
2002 if (present(jobs_opt)) jobs_opt = 0.5d0*jobs_opt
2003
2004 end subroutine oopc_CHobsoperators
2005
2006 !--------------------------------------------------------------------------
2007 ! oopc_obsoperInit
2008 !--------------------------------------------------------------------------
2009 subroutine oopc_obsoperInit(obsSpaceData,headerIndex,columnTrl, &
2010 nmodlev,nobslev,kmode,varno,stnid)
2011 !
2012 !:Purpose: To initialize struct_oopc_obsoperators variables and to allocate
2013 ! arrays.
2014 !
2015 !:Comments:
2016 ! - Allocation of arrays that are dependent on only nlev_bkgrd
2017 ! (nmodlev) have been moved outside this subroutine so that they
2018 ! are allocated only once.
2019 !
2020 !:Arguments:
2021 !
2022 ! :columnTrl: Column of x_background interpolated to observation
2023 ! location. Can have the same vertical levels as the
2024 ! trial field (columnTrl) or as the increment field
2025 ! (columnAnlInc)
2026 ! :kmode: Mode of observation operator
2027 ! - 0 for non-linear/linear model in assimilation (all models
2028 ! included are currently linear)
2029 ! - 1 for determination of sqrt(diag(H*B*H^T))
2030 ! - 2 for tangent linear model
2031 ! - 3 for adjoint model
2032 !
2033 implicit none
2034
2035 ! Arguments:
2036 type(struct_obs), intent(inout) :: obsSpaceData ! Obs-Space Data object
2037 integer, intent(in) :: headerIndex ! Measurement index in obsSpaceData
2038 type(struct_columnData), intent(in) :: columnTrl
2039 integer, intent(in) :: nmodlev ! Number of background field (model) levels
2040 integer, intent(in) :: nobslev ! Number of obs elements (see oopc_obsoper_proceed)
2041 integer, intent(in) :: kmode ! Mode of observation operator
2042 integer, intent(in) :: varno ! obs unit BUFR number
2043 character(len=12), intent(in) :: stnid ! Station ID
2044
2045 ! Locals:
2046 integer :: bodyIndex ! Measurement element index in obsSpaceDate (see oopc_obsoper_proceed)
2047 integer :: jl,nmodlev_uv
2048 real(8), pointer :: col_height_ptr(:)
2049 real(8), allocatable :: uu(:),vv(:)
2050 character(len=2), parameter :: varLevel = 'TH'
2051 real(8) :: checkID
2052
2053 obsoper%nmodlev = nmodlev
2054 obsoper%nobslev = nobslev
2055 obsoper%obs_index = headerIndex
2056 obsoper%varno = varno
2057 obsoper%stnid = stnid
2058
2059 ! Get obs space info that are part of the profile header
2060 obsoper%date = obs_headElem_i(obsSpaceData,OBS_DAT,headerIndex)
2061 obsoper%hhmm = obs_headElem_i(obsSpaceData,OBS_ETM,headerIndex)
2062 ! Constituent identifyer following local version of WMO GRIB Table 08046 (similar to BUFR Table 08043)
2063 obsoper%constituentId = obs_headElem_i(obsSpaceData,OBS_CHM,headerIndex)
2064
2065 ! Check if constituent id is recognized (function will abort if not recognized)
2066 if ( obsoper%constituentId >= 0 .and. obsoper%constituentId < 7000) then
2067 checkID = vnl_varMassFromVarnum(obsoper%constituentId)
2068 end if
2069
2070 obsoper%lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
2071 obsoper%lon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
2072
2073 ! Body info that we only need for first point in the profile
2074 bodyIndex = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2075 ! Index of vertical coordinate type
2076 obsoper%vco = obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex)
2077 ! Model field name (NOMVAR value)
2078 obsoper%varName = vnl_varnameFromVarnum(obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex),obsoper%constituentId,modelName)
2079
2080 ! Allocate arrays
2081
2082 allocate(obsoper%obslev(nobslev)) ! Reference vertical levels
2083 allocate(obsoper%vlayertop(nobslev)) ! Layer tops for layer measurements
2084 allocate(obsoper%vlayerbottom(nobslev)) ! Layer bottoms for layer measurements
2085 allocate(obsoper%modlevindexTop(nobslev)) ! Index of highest model level (lowest index) involved with obs element
2086 allocate(obsoper%modlevindexBot(nobslev)) ! Index of lowest model level (highest index) involved with obs element
2087 allocate(obsoper%obsSpaceTrial(nobslev)) ! Obs estimate from trial field
2088
2089 allocate(obsoper%zh(nobslev,nmodlev)) ! Local model operator H (excluding conversion constants and horizontal interpolation)
2090 allocate(obsoper%zhp(nobslev,nmodlev)) ! Part of zh that excludes aspects related to vertical resolution
2091
2092 obsoper%zh(:,:)=0.0D0
2093 obsoper%zhp(:,:)=0.0D0
2094 obsoper%modlevindexTop(:)=1
2095 obsoper%modlevindexBot(:)=nmodlev
2096
2097 if (.not.col_varExist(columnTrl,'TT')) then
2098 if (oopc_required_field('TT',obsoper%varno)) then
2099 call utl_abort("chem_opsoper_init: TT required for BUFR code " // trim(utl_str(obsoper%varno)))
2100 end if
2101 end if
2102
2103 ! Get background profiles at observation location
2104 do jl=1,nmodlev
2105 obsoper%pp(jl) = col_getPressure(columnTrl,jl,headerIndex,varLevel)
2106 obsoper%trial(jl) = col_getElem(columnTrl,jl,headerIndex,obsoper%varName)
2107 obsoper%tt(jl) = col_getElem(columnTrl,jl,headerIndex,'TT')
2108 end do
2109
2110 if (col_varExist(columnTrl,'TT').and.col_varExist(columnTrl,'HU') &
2111 .and.col_varExist(columnTrl,'P0')) then
2112
2113 ! Height would have been generated in the call to sugomobs.
2114 ! Convert from geopotential to geopotential height.
2115 col_height_ptr => col_getColumn(columnTrl,headerIndex,'Z_T')
2116 obsoper%height(1:nmodlev) = col_height_ptr(1:nmodlev)
2117 else
2118 obsoper%height(:) = -1.
2119 end if
2120
2121 ! Get specific humidity if available
2122 if (col_varExist(columnTrl,'HU')) then
2123 do jl=1,nmodlev
2124 obsoper%hu(jl) = col_getElem(columnTrl,jl,headerIndex,'HU') ! lnq was replaced by q
2125 end do
2126 else
2127 obsoper%hu(:)=-1
2128 end if
2129
2130 ! If applicable, get column upper boundaries for use with total column
2131 ! measurements when the related increment profile is to be restricted to
2132 ! the lower atmosphere (e.g. troposphere or PBL; when tropo_bound>0 )
2133 if (obsoper%vco == 4 .and. nobslev == 1 .and. kmode /= 1) then
2134 if (kmode == 0) then
2135
2136 if (col_varExist(columnTrl,'HU')) then
2137 if (col_varExist(columnTrl,'UU').and.col_varExist(columnTrl,'VV')) then
2138 nmodlev_uv=col_getNumLev(columnTrl,'MM')
2139 allocate(uu(nmodlev_uv),vv(nmodlev_uv))
2140 do jl=1,nmodlev_uv
2141 uu(jl) = col_getElem(columnTrl,jl,headerIndex,'UU')
2142 vv(jl) = col_getElem(columnTrl,jl,headerIndex,'VV')
2143 end do
2144 obsoper%columnBound = oopc_getColBoundary(obsoper%constituentId,&
2145 nmodlev,obsoper%pp,obsoper%tt,obsoper%height, &
2146 hu_opt=obsoper%hu,uu_opt=uu,vv_opt=vv)
2147 deallocate(uu,vv)
2148 else
2149 obsoper%columnBound = oopc_getColBoundary(obsoper%constituentId,nmodlev,obsoper%pp,obsoper%tt,obsoper%height,hu_opt=obsoper%hu)
2150 end if
2151
2152 else
2153 obsoper%columnBound = oopc_getColBoundary(obsoper%constituentId,nmodlev,obsoper%pp,obsoper%tt,obsoper%height)
2154 end if
2155
2156 call oopc_addColBoundary(headerIndex,obsoper%columnBound) ! save boundary for kmode>0 calls using headerIndex
2157
2158 else
2159 obsoper%columnBound = oopc_retrieveColBoundary(headerIndex)
2160 end if
2161 else
2162 obsoper%columnBound = -1.
2163 end if
2164
2165 end subroutine oopc_obsoperInit
2166
2167 !--------------------------------------------------------------------------
2168 ! oopc_obsoper_dealloc
2169 !--------------------------------------------------------------------------
2170 subroutine oopc_obsoper_dealloc
2171 !
2172 !:Purpose: To deallocate arrays for struct_oopc_obsoperators.
2173 !
2174 !:Comments:
2175 ! - Deallocation of arrays that are dependent on only nmodlev have
2176 ! been moved outside this subroutine so that they are
2177 ! deallocated only once.
2178 !
2179 implicit none
2180
2181 if (allocated(obsoper%obslev)) deallocate(obsoper%obslev)
2182 if (allocated(obsoper%vlayertop)) deallocate(obsoper%vlayertop)
2183 if (allocated(obsoper%vlayerbottom)) deallocate(obsoper%vlayerbottom)
2184 if (allocated(obsoper%zh)) deallocate(obsoper%zh)
2185 if (allocated(obsoper%zhp)) deallocate(obsoper%zhp)
2186 if (allocated(obsoper%modlevindexTop)) deallocate(obsoper%modlevindexTop)
2187 if (allocated(obsoper%modlevindexBot)) deallocate(obsoper%modlevindexBot)
2188 if (allocated(obsoper%obsSpaceTrial)) deallocate(obsoper%obsSpaceTrial)
2189 if (allocated(obsoper%success)) deallocate(obsoper%success)
2190 if (allocated(obsoper%ixtr)) deallocate(obsoper%ixtr)
2191
2192 end subroutine oopc_obsoper_dealloc
2193
2194 !---------------------- Routines for observation operators ----------------
2195 !--------------------------------------------------------------------------
2196
2197 !--------------------------------------------------------------------------
2198 ! oopc_obsoperators
2199 !--------------------------------------------------------------------------
2200 subroutine oopc_obsoperators(model_col,obs_col,maxnumHeaders,headerCount,kmode)
2201 !
2202 !:Purpose: To apply observation operator for indicated observation data and
2203 ! condition.
2204 !
2205 !:Arguments:
2206 ! :kmode: mode of observation operator
2207 ! 0) general (potentially non-linear) simulation operator
2208 ! 1) determination of sqrt(diag(H*B*H^T))
2209 ! 2) tangent linear operator
2210 ! 3) linear adjoint operator
2211 !
2212 ! :obs_col: Observation space input/output profile
2213 !
2214 ! +------+------+------+------+------+------+------+
2215 ! |kmode | input/output | profile |
2216 ! +======+====================+====================+
2217 ! | 0 | out | H(xb) |
2218 ! +------+--------------------+--------------------+
2219 ! | 1 | out | sqrt(diag(H*B*H^T))|
2220 ! +------+--------------------+--------------------+
2221 ! | 2 | out | H*dx |
2222 ! +------+--------------------+--------------------+
2223 ! | 3 | in | R**-1 (Hdx-d) |
2224 ! +------+--------------------+--------------------+
2225 !
2226 ! :model_col: Model space input/output profile
2227 !
2228 ! +------+------+------+------+------+------+------+
2229 ! |kmode | input/output | profile |
2230 ! +======+====================+====================+
2231 ! | 0 | in | xb |
2232 ! +------+--------------------+--------------------+
2233 ! | 1 | not used | not used |
2234 ! +------+--------------------+--------------------+
2235 ! | 2 | in | dx at obs location |
2236 ! +------+--------------------+--------------------+
2237 ! | 3 | out | adjoint product |
2238 ! | | | H^T(...) |
2239 ! +------+--------------------+--------------------+
2240 !
2241 ! :maxnumHeaders: Total number of CH obs for this CPU
2242 ! :headerCount: Obs number up to maxnumHeaders
2243 !
2244 ! Further changes required for generalization:
2245 !
2246 ! 1) Add layer average operators.
2247 ! 2) Add AOD operators (summation over model layers).
2248 ! 3) Add option to include use of obs error correlation matrix for kmode=2,3
2249 ! (This may/will need to be done in oop_Hchm and oop_HTchm where the division
2250 ! by stddev_obs is applied. A new routine will be needed for this
2251 ! operation - and others for reading the correlation matrices similarly to the
2252 ! averaging kernels.)
2253 !
2254 ! Comments:
2255 ! - When kmode=0, call from oopc_CHobsoperators passes model_col as a pointer to
2256 ! obsoper%trial.
2257 ! - Does not yet account for potential future applications of obs
2258 ! vertical correlation matrices.
2259 ! - Potential specification of background error std. dev. (fdeStddev(:,1:2)) and correlation matrices
2260 ! for the ensemble-based and lam cases to/could be done when stats for these become in use with constituents.
2261 !
2262 implicit none
2263
2264 ! Arguments:
2265 integer, intent(in) :: kmode
2266 integer, intent(in) :: maxnumHeaders
2267 integer, intent(in) :: headerCount
2268 real(8), intent(inout) :: model_col(obsoper%nmodlev)
2269 real(8), intent(inout) :: obs_col(obsoper%nobslev)
2270
2271 ! Locals:
2272 logical :: successLocal(obsoper%nobslev)
2273 real(8) :: zwork(obsoper%nmodlev),unit_conversion(obsoper%nmodlev)
2274 real(8) :: fdeStddev(obsoper%nmodlev,2),temp_eff(1)
2275 real(8), allocatable :: avg_kern(:,:)
2276 integer :: obslevIndex,modlevIndex,varIndex,nobslevOriginal,varnoOriginal
2277 integer, parameter :: code_len=90
2278 character(len=code_len) :: code ! Must be at least as large as oopc_code_len
2279
2280 if (code_len < oss_obsdata_code_len()) then
2281 call utl_abort('oopc_obsoperators: Length of code string' // &
2282 ' needs to be increased to ' // &
2283 trim(utl_str(oss_obsdata_code_len())))
2284 end if
2285
2286 ! Determine if layer boundaries are assigned to this data source OR, for obsoper%nobslev = 1,
2287 ! if the number of observation operator levels differ (i.e. use of 1D averaging kernels)
2288 ! If so, obtain them for use in this routine.
2289 ! Routine provides obsoper%layerIdentified,
2290 ! obsoper%vlayertop(nobslev),
2291 ! obsoper%vlayerbottom(nobslev)
2292 ! and resets obsoper%nobslev, obsoper%obslev for alternative work vertical levels.
2293 ! The original obsoper%nobslev is saved as nobslevOriginal
2294
2295 nobslevOriginal = obsoper%nobslev
2296 varnoOriginal = obsoper%varno
2297 call oopc_getLevels
2298
2299 ! Prepare observation operator
2300
2301 call oopc_prepareOperator(model_col,avg_kern,nobslevOriginal,maxnumHeaders, &
2302 headerCount,kmode,unit_conversion,successLocal)
2303
2304 ! Finalize required quantities depending on kmode
2305
2306 select case(kmode)
2307
2308 case(0)
2309
2310 ! Finalize non-linear/linear operator step
2311
2312 if (obsoper%iavgkern > 0) then
2313 if (oopc_checkType(oopc_avgkern%stnids, &
2314 oopc_avgkern%type,obsoper%stnid,'log')) then
2315
2316 zwork(1:obsoper%nobslev) = 0.0d0
2317 do obslevIndex=1,obsoper%nobslev
2318 if (successLocal(obslevIndex)) &
2319 zwork(obslevIndex)=log(dot_product(obsoper%zh(obslevIndex, &
2320 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)), &
2321 model_col(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) ))
2322 end do
2323
2324 do obslevIndex=1,oopc_avgkern%n_lvl(obsoper%iavgkern)
2325 if (obsoper%success(obslevIndex)) then
2326 obs_col(obslevIndex)= dot_product(avg_kern(obslevIndex, &
2327 1:obsoper%nobslev),zwork(1:obsoper%nobslev))
2328 end if
2329
2330 ! Add a priori contribution when provided
2331 if (oopc_avgkern%n_col(obsoper%iavgkern) >= obsoper%nobslev+1) then
2332 obs_col(obslevIndex) = obs_col(obslevIndex) + &
2333 avg_kern(obslevIndex,obsoper%nobslev+1)
2334 end if
2335
2336 obs_col(obslevIndex) = exp(obs_col(obslevIndex))
2337 end do
2338
2339 if ( oopc_avgkern%n_col(obsoper%iavgkern) == obsoper%nobslev+2 &
2340 .and. nobslevOriginal == 1 ) then
2341
2342 obs_col(1) = sum( avg_kern(1:oopc_avgkern%n_lvl(obsoper%iavgkern), &
2343 obsoper%nobslev+2)*obs_col(1:oopc_avgkern%n_lvl(obsoper%iavgkern)) )
2344 end if
2345 else
2346
2347 ! Standard treatment for linear model
2348
2349 do obslevIndex=1,oopc_avgkern%n_lvl(obsoper%iavgkern)
2350 if (obsoper%success(obslevIndex)) then
2351 obs_col(obslevIndex)=dot_product(obsoper%zh(obslevIndex,&
2352 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)), &
2353 model_col(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)))
2354 end if
2355 ! Add a priori contribution when provided
2356 if ( oopc_avgkern%n_col(obsoper%iavgkern) >= obsoper%nobslev+1) then
2357 obs_col(obslevIndex) = obs_col(obslevIndex) + &
2358 avg_kern(obslevIndex,obsoper%nobslev+1)
2359 end if
2360 end do
2361
2362 ! Account for integration via weighted summation
2363 if ( oopc_avgkern%n_col(obsoper%iavgkern) == obsoper%nobslev+2 &
2364 .and. nobslevOriginal == 1 .and. &
2365 oopc_avgkern%n_lvl(obsoper%iavgkern) > 1 ) then
2366
2367 obs_col(1) = sum( avg_kern(1:oopc_avgkern%n_lvl(obsoper%iavgkern), &
2368 obsoper%nobslev+2)*obs_col(1:oopc_avgkern%n_lvl(obsoper%iavgkern)) )
2369 end if
2370 end if
2371
2372 else
2373
2374 ! Standard treatment for linear model
2375
2376 do obslevIndex=1,nobslevOriginal
2377 if (obsoper%success(obslevIndex)) then
2378 obs_col(obslevIndex)=dot_product(obsoper%zh(obslevIndex, &
2379 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)), &
2380 model_col(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)))
2381 end if
2382 end do
2383
2384 end if
2385
2386 ! Calculate concentration-weighted effective temperature (for output purpose)
2387
2388 if ((obsoper%constituentId >= 0 .and. obsoper%constituentId < 7000) .and. &
2389 trim(obsoper%operatorCategory) == 'Integ' .and. obsoper%nobslev == 1 .and. &
2390 obsoper%vco == 4 .and. obsoper%success(1)) then
2391
2392 if (all(obsoper%tt > 0.0) .and. obs_col(1) > 0.0) then
2393 temp_eff(1)=dot_product(obsoper%zh(1, &
2394 obsoper%modlevindexTop(1):obsoper%modlevindexBot(1)), &
2395 obsoper%tt(obsoper%modlevindexTop(1):obsoper%modlevindexBot(1))* &
2396 model_col(obsoper%modlevindexTop(1):obsoper%modlevindexBot(1))) &
2397 /obs_col(1)
2398 code=oss_obsdata_get_header_code(obsoper%lon,obsoper%lat,&
2399 obsoper%date,obsoper%hhmm,obsoper%stnid)
2400 call oopc_addEfftempObsdata(code,temp_eff)
2401 end if
2402 end if
2403
2404 case(1)
2405
2406 ! Compute sqrt(diag(H*B*H^T))
2407
2408 ! Apply unit conversion to observation operator
2409 do obslevIndex=1,nobslevOriginal
2410 if (obsoper%success(obslevIndex)) then
2411 obsoper%zh(obslevIndex,:) = unit_conversion * obsoper%zh(obslevIndex,:)
2412 end if
2413 end do
2414
2415 ! Get background error std dev profile at obs location
2416 fdeStddev(:,:)=0
2417 call bcsc_getBgStddev(obsoper%varName,obsoper%nmodlev, &
2418 obsoper%lat,obsoper%lon,fdeStddev(:,1))
2419
2420 ! Identify variable position index in background error correlation matrices
2421
2422 varIndex=1
2423 do while (trim(bgStats%varNameList(varIndex)) /= '')
2424 if (trim(bgStats%varNameList(varIndex)) == trim(obsoper%varName)) exit
2425 varIndex=varIndex+1
2426 end do
2427 if (trim(bgStats%varNameList(varIndex)) == '') then
2428 call utl_abort('oopc_obsoperators: Correlation matrix not found for ' &
2429 // trim(obsoper%varName) )
2430 end if
2431
2432 do obslevIndex=1,nobslevOriginal
2433 if (obsoper%success(obslevIndex)) then
2434 do modlevIndex=obsoper%modlevindexTop(obslevIndex), &
2435 obsoper%modlevindexBot(obslevIndex)
2436
2437 zwork(modlevIndex)=sum(obsoper%zh(obslevIndex, &
2438 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) &
2439 *bgStats%corvert(modlevIndex,obsoper%modlevindexTop(obslevIndex): &
2440 obsoper%modlevindexBot(obslevIndex),varIndex) &
2441 *fdeStddev(obsoper%modlevindexTop(obslevIndex): &
2442 obsoper%modlevindexBot(obslevIndex),1)) &
2443 *fdeStddev(modlevIndex,1)
2444 end do
2445 obs_col(obslevIndex)=sum(obsoper%zh(obslevIndex, &
2446 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) &
2447 *zwork(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)))
2448
2449 obs_col(obslevIndex) = sqrt(obs_col(obslevIndex)) ! save as sqrt(h*B*h^T)
2450 else
2451 obs_col(obslevIndex) = 0.0
2452 end if
2453 end do
2454
2455 case(2)
2456
2457 ! Finalize TL operator application
2458
2459 ! Standard treatment for linear and linearized model
2460
2461 if (obsoper%iavgkern > 0) then
2462
2463 do obslevIndex=1,oopc_avgkern%n_lvl(obsoper%iavgkern)
2464 if (obsoper%success(obslevIndex)) then
2465 obs_col(obslevIndex)=dot_product(obsoper%zh(obslevIndex, &
2466 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)), &
2467 model_col(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)))
2468 end if
2469 end do
2470
2471 ! Account for integration via weighted summation
2472 if ( oopc_avgkern%n_col(obsoper%iavgkern) == obsoper%nobslev+2 &
2473 .and. nobslevOriginal == 1 .and. &
2474 oopc_avgkern%n_lvl(obsoper%iavgkern) > 1 ) then
2475
2476 obs_col(1) = sum( avg_kern(1:oopc_avgkern%n_lvl(obsoper%iavgkern), &
2477 obsoper%nobslev+2)*obs_col(1:oopc_avgkern%n_lvl(obsoper%iavgkern)) )
2478
2479 end if
2480 else
2481
2482 do obslevIndex=1,nobslevOriginal
2483 if (obsoper%success(obslevIndex)) then
2484 obs_col(obslevIndex)=dot_product(obsoper%zh(obslevIndex, &
2485 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)), &
2486 model_col(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)))
2487 else
2488 obs_col(obslevIndex)=0.0d0
2489 end if
2490 end do
2491
2492 end if
2493
2494 case(3)
2495
2496 ! H^T*grad contribution from adjoint of tangent linear model.
2497
2498 model_col(:) = 0.0
2499
2500 ! Standard treatment for linear and linearized model
2501
2502 if (obsoper%iavgkern > 0) then
2503
2504 if ( oopc_avgkern%n_col(obsoper%iavgkern) == obsoper%nobslev+2 &
2505 .and. nobslevOriginal == 1 .and. &
2506 oopc_avgkern%n_lvl(obsoper%iavgkern) > 1 ) then
2507
2508 ! Account for integration via weighted summation
2509 zwork(:)=0.0d0
2510 do obslevIndex=obsoper%modlevindexTop(obslevIndex),obsoper%modlevindexBot(obslevIndex)
2511 zwork(obslevIndex)= &
2512 sum(obsoper%zh(1:oopc_avgkern%n_lvl(obsoper%iavgkern),obslevIndex)* &
2513 avg_kern(1:oopc_avgkern%n_lvl(obsoper%iavgkern),obsoper%nobslev+2) )
2514 end do
2515 obsoper%zh(1,:) = 0.0d0
2516 obsoper%zh(1,obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) = &
2517 zwork(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex))
2518 end if
2519 end if
2520
2521 do obslevIndex=1,nobslevOriginal
2522 if (obsoper%success(obslevIndex)) then
2523 zwork(:)=0.0D0
2524 zwork(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) = &
2525 obs_col(obslevIndex)*obsoper%zh(obslevIndex, &
2526 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex))
2527
2528 call oopc_convertUnits(zwork,incr_opt=.true.)
2529
2530 model_col(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) = &
2531 model_col(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) + &
2532 zwork(obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex))
2533 end if
2534 end do
2535
2536 end select
2537
2538 if (obsoper%iavgkern > 0) deallocate(avg_kern)
2539
2540 end subroutine oopc_obsoperators
2541
2542 !--------------------------------------------------------------------------
2543 ! oopc_prepareOperator
2544 !--------------------------------------------------------------------------
2545 subroutine oopc_prepareOperator(model_col,avg_kern,nobslevOriginal, &
2546 maxnumHeaders,headerCount,kmode,unit_conversion,successLocal)
2547 !
2548 !:Purpose: To prepare observation operator
2549 !
2550 implicit none
2551
2552 ! Arguments:
2553 integer, intent(in) :: kmode ! Mode of observation operator
2554 integer, intent(in) :: nobslevOriginal ! Number of actual obs elements in measurement
2555 integer, intent(in) :: maxnumHeaders ! Total number of CH obs for this CPU
2556 integer, intent(in) :: headerCount ! Obs counter/index
2557 real(8), intent(inout) :: model_col(obsoper%nmodlev)
2558 real(8), intent(out) :: unit_conversion(obsoper%nmodlev)
2559 logical, intent(out) :: successLocal(obsoper%nobslev)
2560 real(8), allocatable, intent(out) :: avg_kern(:,:)
2561
2562 ! Locals:
2563 real(8), allocatable :: press_obs(:)
2564 integer, allocatable :: ixtrLocal(:)
2565 real(8) :: fdeStddev(obsoper%nmodlev,2),zhmin
2566 integer :: obslevIndex,modlevIndex
2567 integer, parameter :: code_len=90
2568 character(len=code_len) :: code ! Must be at least as large as oopc_code_len
2569 character(len=20) :: message
2570 logical :: successDepot
2571
2572 ! Check if obs BUFR element needs to be changed with use of averaging kernels
2573
2574 if ( obsoper%iavgkern > 0 ) then
2575 if ( oopc_avgkern%n_col(obsoper%iavgkern) == obsoper%nobslev+2 &
2576 .and. nobslevOriginal == 1 ) then
2577
2578 obsoper%varno=oopc_avgkern%ProfElement(obsoper%iavgkern)
2579 end if
2580 end if
2581
2582 ! Identify observation operator based on observation units and presence or
2583 ! not of layer boundaries
2584
2585 if (bufr_IsIntegral(obsoper%varno)) then
2586 if (.not.obsoper%layerIdentified) then
2587 write(*,*) '----------------------------------------------------------'
2588 write(*,*) 'STNID, BUFR index, nobslev: ',obsoper%stnid,' ', &
2589 obsoper%varno,obsoper%nobslev
2590 call utl_abort('oopc_obsoperators: Required layer boundaries not available!')
2591 else
2592 ! Vertical integration operator
2593 obsoper%operatorCategory='Integ'
2594 end if
2595 else if (obsoper%layerIdentified) then
2596 ! Layer averaging operator
2597 obsoper%operatorCategory='LayerAvg'
2598 else if (obsoper%vco == 5 .and. obsoper%nobslev == 1) then
2599 ! Surface point (in-situ) measurement
2600 obsoper%operatorCategory='Surface'
2601 else
2602 ! Vertical interpolation operator
2603 obsoper%operatorCategory='Interp'
2604 end if
2605
2606 ! Look for pre-calculated operator
2607
2608 successDepot = oopc_operatorDepot(headerCount,maxnumHeaders,kmode,'get')
2609
2610 if ( successDepot ) then
2611
2612 ! Get averaging kernels if requested
2613
2614 if (obsoper%iavgkern > 0) then
2615
2616 if (allocated(avg_kern)) deallocate(avg_kern)
2617 allocate(avg_kern(oopc_avgkern%n_lvl(obsoper%iavgkern), &
2618 oopc_avgkern%n_col(obsoper%iavgkern)))
2619
2620 code=oss_obsdata_get_header_code(obsoper%lon,obsoper%lat,obsoper%date, &
2621 obsoper%hhmm,obsoper%stnid)
2622 call oopc_getAvgkern(obsoper%iavgkern,oopc_avgkern%n_lvl(obsoper%iavgkern), &
2623 oopc_avgkern%n_col(obsoper%iavgkern),code,avg_kern)
2624 end if
2625
2626 ! Apply unit conversion (apply later when kmode=3)
2627 if (kmode == 2) call oopc_convertUnits(model_col,incr_opt=.true.)
2628
2629 return
2630 end if
2631
2632 ! Indicate if the generalized innovation operator is to be applied.
2633
2634 obsoper%applyGenOper=.false.
2635 if (obsoper%constituentId >= 0 .and. &
2636 oopc_checkType(operatorSubType(1,:),operatorSubType(2,:), &
2637 obsoper%stnid,'genOper')) then
2638
2639 if (kmode /= 1 .and. (trim(obsoper%operatorCategory) == 'Integ' .or. &
2640 trim(obsoper%operatorCategory) == 'LayerAvg' )) then
2641
2642 if ( kmode == 2) then
2643 ! Set reference profiles for use with generalized innovation operator
2644 ! when kmode>=2
2645 call oopc_addToProfileSet(oopc_climFields,oopc_bgRef, &
2646 oopc_constituentsSize,2,obsoper%nmodlev,obsoper%pp, &
2647 obsoper%height,obsoper%lat,obsoper%lon, obsoper%obs_index, &
2648 oopc_obsdata_maxsize,varKind_opt='CH', &
2649 varNumber_opt=obsoper%constituentId, &
2650 tt_opt=obsoper%tt,hu_opt=obsoper%hu)
2651
2652 ! Get background error std dev profile at obs locations
2653 fdeStddev(:,:)=0.D0
2654
2655 call bcsc_getBgStddev(obsoper%varName,obsoper%nmodlev,obsoper%lat, &
2656 obsoper%lon,fdeStddev(:,1),vlev_opt=obsoper%pp)
2657
2658 call bcsc_addBgStddev(obsoper%obs_index,fdeStddev, &
2659 oopc_obsdata_maxsize)
2660
2661 end if
2662 if (kmode >= 2) obsoper%applyGenOper = .true.
2663
2664 end if
2665 end if
2666
2667 ! Apply unit conversion (apply unit conversion later for kmode=3)
2668
2669 select case(kmode)
2670 case(0)
2671 ! Perform transformation and unit conversion on obsoper%trial.
2672 ! Note that model_col => obsoper%trial for kmode=0.
2673 call oopc_convertUnits(obsoper%trial)
2674 case(1)
2675 ! Save the conversion factor in <unit_conversion>
2676 unit_conversion(:) = 1.0
2677 call oopc_convertUnits(unit_conversion,incr_opt=.true.)
2678 case(2)
2679 call oopc_convertUnits(model_col,incr_opt=.true.)
2680 end select
2681
2682 if (obsoper%applyGenOper) then
2683 ! Perform unit conversion on obsoper%trial when applying the generalized
2684 ! obs operator for kmode=2,3. Keep obsoper%trial in ug/kg in this case.
2685 call oopc_convertUnits(obsoper%trial,ppb_opt=.true.)
2686 end if
2687
2688 if (allocated(press_obs)) deallocate(press_obs,ixtrLocal)
2689 allocate(press_obs(obsoper%nobslev),ixtrLocal(obsoper%nobslev))
2690 if (nobslevOriginal == obsoper%nobslev) then
2691 ixtrLocal(:) = obsoper%ixtr(:)
2692 successLocal(:) = obsoper%success(:)
2693 else
2694 if (any(obsoper%ixtr(:) == 0)) then
2695 ixtrLocal(:)=0
2696 else
2697 ixtrLocal(:)=1
2698 end if
2699 if (any(obsoper%success)) then
2700 successLocal(:) = .true.
2701 else
2702 successLocal(:) = .false.
2703 end if
2704 end if
2705
2706 ! Convert observation vertical coordinate value(s) to pressure if needed
2707
2708 select case(obsoper%vco)
2709 case(1)
2710 ! Convert altitude to pressure
2711 if (trim(obsoper%operatorCategory) == 'Interp') then
2712 press_obs = phf_convert_z_to_pressure(obsoper%obslev,obsoper%height, &
2713 obsoper%pp, &
2714 obsoper%nobslev,obsoper%nmodlev,obsoper%lat,successLocal)
2715 ! Allows for obs levels below the lowest TH level and above the surface
2716 where (obsoper%obslev(1:obsoper%nobslev) < &
2717 obsoper%height(obsoper%nmodlev)) &
2718 press_obs(1:obsoper%nobslev)= obsoper%pp(obsoper%nmodlev)
2719
2720 else if (trim(obsoper%operatorCategory) == 'Integ' .or. &
2721 trim(obsoper%operatorCategory) == 'LayerAvg') then
2722 obsoper%vlayertop = phf_convert_z_to_pressure(obsoper%vlayertop, &
2723 obsoper%height, &
2724 obsoper%pp,obsoper%nobslev,obsoper%nmodlev, &
2725 obsoper%lat,successLocal)
2726 obsoper%vlayerbottom = phf_convert_z_to_pressure(obsoper%vlayerbottom, &
2727 obsoper%height, &
2728 obsoper%pp,obsoper%nobslev,obsoper%nmodlev, &
2729 obsoper%lat,successLocal)
2730 end if
2731 case(2)
2732 ! Pressure, no conversion needed
2733 if (trim(obsoper%operatorCategory) == 'Interp') press_obs = obsoper%obslev
2734 case(4,5)
2735 ! No actions taken
2736 case default
2737 call utl_abort("oopc_obsoperators: vertical coordinate type vco = " &
2738 // trim(utl_str(obsoper%vco)) // &
2739 " not available for this operator.")
2740 end select
2741
2742 ! Determine if averaging kernel is to be applied
2743
2744 if (nobslevOriginal == obsoper%nobslev) then
2745 obsoper%iavgkern = oopc_findAvgkern(obsoper%stnid,obsoper%varno, &
2746 obsoper%nobslev)
2747 else
2748 obsoper%iavgkern = oopc_findAvgkern(obsoper%stnid,obsoper%varno,0)
2749 end if
2750
2751 ! Apply appropriate core observation operator
2752
2753 select case(trim(obsoper%operatorCategory))
2754 case('Interp')
2755
2756 ! Vertical interpolation operator
2757
2758 if ( obsoper%iavgkern /= 0 ) then
2759 message = 'doAll&noExtrap'
2760 else
2761 message = 'noExtrap'
2762 end if
2763 call ppo_vertInterpWgts(obsoper%pp,press_obs,obsoper%nmodlev, &
2764 obsoper%nobslev,obsoper%zh,obsoper%modlevindexTop, &
2765 obsoper%modlevindexBot,method_opt=oopc_getType(operatorSubType(1,:), &
2766 operatorSubType(2,:),obsoper%stnid),skipType_opt=message, &
2767 outbound_opt=ixtrLocal,success_opt=successLocal)
2768
2769 case('Surface')
2770
2771 ! Surface point measurement
2772
2773 ! Set weight to unity for lowest model level.
2774 obsoper%zh(1,1:obsoper%nmodlev-1)=0.0
2775 obsoper%zh(1,obsoper%nmodlev) = 1.0
2776
2777 ! Set range of elements for model vertical levels
2778 obsoper%modlevindexTop(1) = obsoper%nmodlev
2779 obsoper%modlevindexBot(1) = obsoper%nmodlev
2780
2781 case('Integ')
2782
2783 ! Layer integration operator
2784
2785 if ( obsoper%iavgkern /= 0 ) then
2786 message = 'doAll&noExtrap'
2787 else
2788 message = 'default'
2789 end if
2790
2791 call oopc_vertObsLayersWgts('integ',ixtrLocal,successLocal,kmode,message)
2792
2793 case('LayerAvg')
2794
2795 ! Layer averaging operator
2796
2797 if ( obsoper%iavgkern /= 0 ) then
2798 message = 'doAll&noExtrap'
2799 else
2800 message = 'default'
2801 end if
2802
2803 call oopc_vertObsLayersWgts('avg',ixtrLocal,successLocal,kmode,message)
2804
2805 end select
2806
2807 if (nobslevOriginal == obsoper%nobslev) then
2808 obsoper%ixtr(:) = ixtrLocal(:)
2809 obsoper%success(:) = successLocal(:)
2810 else
2811 if (all(ixtrLocal(:) == 0)) then
2812 obsoper%ixtr(:)=0
2813 else
2814 obsoper%ixtr(:)=1
2815 end if
2816 if (all(successLocal)) then
2817 obsoper%success(:) = .true.
2818 else
2819 obsoper%success(:) = .false.
2820 end if
2821 end if
2822 deallocate(press_obs,ixtrLocal)
2823
2824 ! Apply averaging kernels if requested
2825
2826 if (obsoper%iavgkern > 0) then
2827
2828 allocate(avg_kern(oopc_avgkern%n_lvl(obsoper%iavgkern), &
2829 oopc_avgkern%n_col(obsoper%iavgkern)))
2830
2831 code=oss_obsdata_get_header_code(obsoper%lon,obsoper%lat, &
2832 obsoper%date,obsoper%hhmm,obsoper%stnid)
2833 call oopc_getAvgkern(obsoper%iavgkern, &
2834 oopc_avgkern%n_lvl(obsoper%iavgkern), &
2835 oopc_avgkern%n_col(obsoper%iavgkern),code,avg_kern)
2836
2837 if (oopc_checkType(oopc_avgkern%stnids,oopc_avgkern%type,obsoper%stnid, &
2838 'default')) then
2839 do obslevIndex=1,oopc_avgkern%n_lvl(obsoper%iavgkern)
2840 if (obsoper%success(obslevIndex)) then
2841
2842 ! Apply averaging kernels to observation operator(s)
2843
2844 obsoper%zh(obslevIndex,:) = matmul(avg_kern(obslevIndex, &
2845 1:obsoper%nobslev),obsoper%zh(:,:))
2846 if (obsoper%applyGenOper) then
2847 obsoper%zhp(obslevIndex,:) = &
2848 matmul(avg_kern(obslevIndex,1:obsoper%nobslev),obsoper%zhp(:,:))
2849 end if
2850
2851 ! Extend vertical range of obs operator according to the influence of
2852 ! the averaging kernel. Either extend to the entire model vertical range
2853 ! (commented out below) or to the vertical range with non-negligable values.
2854
2855 ! obsoper%modlevindexTop(obslevIndex) = 1
2856 ! obsoper%modlevindexBot(obslevIndex) = obsoper%nmodlev
2857
2858 zhmin=1.0D-10*maxval(abs(obsoper%zh(obslevIndex,:)))
2859 do modlevIndex=1,obsoper%modlevindexTop(obslevIndex)
2860 if (abs(obsoper%zh(obslevIndex,modlevIndex)) > zhmin) exit
2861 end do
2862 if (modlevIndex > obsoper%modlevindexTop(obslevIndex)) then
2863 modlevIndex=obsoper%modlevindexTop(obslevIndex)
2864 end if
2865 obsoper%modlevindexTop(obslevIndex) = modlevIndex
2866 do modlevIndex=obsoper%nmodlev,obsoper%modlevindexBot(obslevIndex),-1
2867 if (abs(obsoper%zh(obslevIndex,modlevIndex)) > zhmin) exit
2868 end do
2869 if (modlevIndex.lt.obsoper%modlevindexBot(obslevIndex)) then
2870 modlevIndex=obsoper%modlevindexBot(obslevIndex)
2871 end if
2872 obsoper%modlevindexBot(obslevIndex) = modlevIndex
2873
2874 end if
2875 end do
2876
2877 else if (oopc_checkType(oopc_avgkern%stnids,oopc_avgkern%type, &
2878 obsoper%stnid,'log')) then
2879
2880 if (kmode == 0) then
2881
2882 ! Apply log-space averaging kernel below - no transformation needed here
2883 ! Do not merge the averaging kernel (avgkern) and the vertical interpolator (zh)
2884
2885 else
2886
2887 ! Apply linearization of operator involving the log-space averaging kernel
2888
2889 do obslevIndex=1,obsoper%nobslev
2890 avg_kern(1:oopc_avgkern%n_lvl(obsoper%iavgkern),obslevIndex) = &
2891 avg_kern(1:oopc_avgkern%n_lvl(obsoper%iavgkern),obslevIndex) / &
2892 dot_product( obsoper%zh(obslevIndex, &
2893 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)), &
2894 obsoper%trial(obsoper%modlevindexTop(obslevIndex): &
2895 obsoper%modlevindexBot(obslevIndex)) )
2896 end do
2897
2898 do obslevIndex=1,oopc_avgkern%n_lvl(obsoper%iavgkern)
2899 if (obsoper%success(obslevIndex)) then
2900 avg_kern(obslevIndex,1:obsoper%nobslev) = &
2901 obsoper%obsSpaceTrial(obslevIndex)*avg_kern(obslevIndex, &
2902 1:obsoper%nobslev)
2903
2904 ! Merge the averaging kernel matrix (avgkern) and the
2905 ! vertical interpolator (initial zh)
2906 obsoper%zh(obslevIndex,:) = matmul(avg_kern(obslevIndex, &
2907 1:obsoper%nobslev),obsoper%zh(:,:))
2908
2909 end if
2910 end do
2911
2912 end if
2913
2914 ! Note that obsoper%applyGenOper=.true. is not set up for this case
2915 if (obsoper%applyGenOper) then
2916 call utl_abort('prepareOperator: Log ' // &
2917 'space averaging kernels not currently usable with ' // &
2918 'obsoper%applyGenOper=.true.')
2919 end if
2920 else
2921 call utl_abort('oopc_prepareOperator: This averaging kernel ' // &
2922 'application not yet available')
2923 end if
2924 else
2925 if ( nobslevOriginal /= obsoper%nobslev ) then
2926 call utl_abort('oops_prepareOperator: Case of differing obs ' // &
2927 'and calculation levels not recognized.')
2928 end if
2929 end if
2930
2931 ! Apply generalized innovation operator if requested
2932
2933 if (obsoper%applyGenOper) call oopc_genOper(kmode)
2934
2935 ! Save operator if needed
2936
2937 successDepot = oopc_operatorDepot(headerCount,maxnumHeaders,kmode,'save')
2938
2939 end subroutine oopc_prepareOperator
2940
2941 !--------------------------------------------------------------------------
2942 ! oopc_operatorDepot
2943 !--------------------------------------------------------------------------
2944 function oopc_operatorDepot(headerCount,maxnumHeaders,kmode,action) result(success)
2945 !
2946 !:Purpose: To save or get previously calculated linear observation operator
2947 ! for use by TL and AD operators
2948 !
2949 implicit none
2950
2951 ! Arguments:
2952 integer, intent(in) :: kmode ! Mode of observation operator
2953 integer, intent(in) :: headerCount ! Obs counter/index
2954 integer, intent(in) :: maxnumHeaders ! Total number of CH obs for this CPU
2955 character(len=*), intent(in) :: action ! 'save' or 'get' operator for current obs
2956 ! Result:
2957 logical :: success ! Succes of action succeeded
2958
2959 ! Locals:
2960 integer, save :: maxnumOperators
2961 type(struct_oopc_operatorsDepot), allocatable, save :: operators(:)
2962
2963 success = .false.
2964
2965 if ( .not.oopc_storeOperators .or. kmode < 2 .or. headerCount <= 0 .or. &
2966 ( any(oopc_tropo_mode(:) >= 1) .and. &
2967 trim(obsoper%operatorCategory) == 'Integ') ) return
2968
2969 if ( trim(action) == 'get' .and. initializedOperators .and. kmode >=2 ) then
2970
2971 ! Check for available pre-calculated operator
2972
2973 if (operators(headerCount)%nobslev > 0) then
2974
2975 ! Update arrays/values that might have changed or would otherwise be set
2976 ! below when not stored earlier.
2977
2978 obsoper%trial(:) = operators(headerCount)%trial(:)
2979 obsoper%vlayertop(:) = operators(headerCount)%vlayertop(:)
2980 obsoper%vlayerbottom(:) = operators(headerCount)%vlayerbottom(:)
2981 obsoper%modlevindexTop(:) = operators(headerCount)%modlevindexTop(:)
2982 obsoper%modlevindexBot(:) = operators(headerCount)%modlevindexBot(:)
2983 obsoper%zh(:,:) = operators(headerCount)%zh(:,:)
2984 obsoper%ixtr(:) = operators(headerCount)%ixtr(:)
2985 obsoper%success(:) = operators(headerCount)%success(:)
2986 obsoper%iavgkern = operators(headerCount)%iavgkern
2987 obsoper%applyGenOper = operators(headerCount)%applyGenOper
2988
2989 if ( obsoper%applyGenOper ) then
2990 obsoper%zhp(:,:) = operators(headerCount)%zhp(:,:)
2991 end if
2992
2993 success = .true.
2994 else
2995 ! Calculations will be performed in th calling routine to prepare
2996 ! the operator fields
2997 end if
2998
2999 else if ( trim(action) == 'save' .and. kmode == 2 ) then
3000
3001 ! Save operator for later use
3002
3003 if ( .not.initializedOperators ) then
3004 if (allocated(operators)) deallocate(operators)
3005 maxnumOperators = maxnumHeaders
3006 allocate(operators(maxnumOperators))
3007 initializedOperators = .true.
3008
3009 operators(:)%nobslev = 0
3010
3011 write(*,*) 'oopc_operatorDepot: Max number of operators to save: ', &
3012 maxnumOperators
3013
3014 end if
3015
3016 if (allocated(operators(headerCount)%zh)) then
3017 deallocate(operators(headerCount)%trial)
3018 deallocate(operators(headerCount)%vlayertop)
3019 deallocate(operators(headerCount)%vlayertop)
3020 deallocate(operators(headerCount)%modlevindexTop)
3021 deallocate(operators(headerCount)%modlevindexBot)
3022 deallocate(operators(headerCount)%zh)
3023 deallocate(operators(headerCount)%ixtr)
3024 deallocate(operators(headerCount)%success)
3025 if ( obsoper%applyGenOper ) deallocate(operators(headerCount)%zhp)
3026 end if
3027
3028 allocate(operators(headerCount)%trial(obsoper%nmodlev))
3029 allocate(operators(headerCount)%vlayertop(obsoper%nobslev))
3030 allocate(operators(headerCount)%vlayerbottom(obsoper%nobslev))
3031 allocate(operators(headerCount)%modlevindexTop(obsoper%nobslev))
3032 allocate(operators(headerCount)%modlevindexBot(obsoper%nobslev))
3033 allocate(operators(headerCount)%zh(obsoper%nobslev,obsoper%nmodlev))
3034 allocate(operators(headerCount)%success(obsoper%nobslev))
3035 allocate(operators(headerCount)%ixtr(obsoper%nobslev))
3036
3037 operators(headerCount)%nobslev=obsoper%nobslev
3038 operators(headerCount)%iavgkern=obsoper%iavgkern
3039 operators(headerCount)%applyGenOper=obsoper%applyGenOper
3040 operators(headerCount)%trial(:) = obsoper%trial(:)
3041 operators(headerCount)%vlayertop(:) = obsoper%vlayertop(:)
3042 operators(headerCount)%vlayerbottom(:) = obsoper%vlayerbottom(:)
3043 operators(headerCount)%modlevindexTop(:) = obsoper%modlevindexTop(:)
3044 operators(headerCount)%modlevindexBot(:) = obsoper%modlevindexBot(:)
3045 operators(headerCount)%zh(:,:) = obsoper%zh(:,:)
3046 operators(headerCount)%ixtr(:) = obsoper%ixtr(:)
3047 operators(headerCount)%success(:) = obsoper%success(:)
3048
3049 if ( obsoper%applyGenOper ) then
3050 allocate(operators(headerCount)%zhp(obsoper%nobslev,obsoper%nmodlev))
3051 operators(headerCount)%zhp(:,:) = obsoper%zhp(:,:)
3052 end if
3053
3054 success = .true.
3055
3056 end if
3057
3058 end function oopc_operatorDepot
3059
3060 !--------------------------------------------------------------------------
3061 ! oopc_convertUnits
3062 !--------------------------------------------------------------------------
3063 subroutine oopc_convertUnits(model_col,ppb_opt,incr_opt)
3064 !
3065 !:Purpose: To set unit-conversion factor for consistency of Hx units with
3066 ! obs units.
3067 !
3068 !:Arguments:
3069 ! :ppb_opt: indicates whether model_col should be kept in
3070 ! ug/kg instead of the units dictated by the BUFR
3071 ! number (optional, .false. by default)
3072 ! :incr_opt: indicates if model_col is actually an increment
3073 ! (optional, .false. by default). Needed for non-linear
3074 ! transformations (i.e. for 'HU')
3075 ! :model_col Array to be converted. Either trial or increment-related.
3076 !
3077 !:Comments:
3078 !
3079 ! A. Standard model/analysis species field provided as mass mixing
3080 ! ratio in ug/kg (ppb). Conversion to ppb is applied when this is
3081 ! not the case except for AOD and surface emissions.
3082 ! As this is hard-coded, any changes in analysis variable must
3083 ! be reflected by correspondingly modifying this module.
3084 !
3085 ! B. Unit-conversion factor is calculated in oopc_convertUnits
3086 ! from the following factors:
3087 ! (1) physical constants
3088 ! (2) parameters related to a particular species such as molecular
3089 ! mass
3090 ! (3) variables such as T and P from background field at each
3091 ! iteration
3092 !
3093 ! C. The baseline integral observation operator can be interpreted as
3094 ! being integrals of the gas partial pressure, giving products in
3095 ! kg/m^2, e.g. with sample discretized layer integrals::
3096 ! (mass density) * dz = - d(gas partial pressure)/g
3097 ! = - [rho(gas)/rho(air)]*dP/g
3098 ! = - 1E-9 * [mass mixing ratio in parts per billion (ppb)]*dP/g
3099 !
3100 ! The actual integration in pressure (in Pascal) is performed
3101 ! outside this routine. For integral products in kg/m^2, the output
3102 ! of this routine is to be in mmr/(m/s^2) (mmr=mass mixing ratio),
3103 ! which is equivalent to (1E-9 ug/kg)/(m/s^2) and kg/(m^2*Pa).
3104 ! Therefore, the input value in ug/kg has to be multiplied by 1E-9/g
3105 ! (g=RG below).
3106 !
3107 ! For integral products in other units, additional conversion
3108 ! factors are also to be applied.
3109 !
3110 ! D. List should be revised following changes to the 'tableburp' file.
3111 !
3112 !
3113 ! E. Coefficients related to unit conversion
3114 !
3115 ! rho_stp=1.293 Air density at STP (1.293 kg/m^3)
3116 ! RG=9.807 (=g) Acceleration due to gravity (m/s^2)
3117 ! MPC_AVOGADRO_R8 = Na Avogadro's number. 6.023E23 molecules/mole
3118 ! MPC_MOLAR_MASS_DRY_AIR_R8 (m_air) Dry air molecular mass. 28.9644 g/mole
3119 ! MPC_RGAS_IDEAL_R8 = R Ideal gas constant. 8.341 J/mole/K (J=kg m^2/s^2)
3120 !
3121 ! PV = nRT (n=number of moles)
3122 !
3123 ! MPC_RGAS_DRY_AIR_R8 = Rd Dry air constant. 287.1 J/kg/K (J=kg m^2/s^2)
3124 ! = MPC_RGAS_IDEAL_R8 * 1000 g/km / MPC_MOLAR_MASS_DRY_AIR_R8
3125 !
3126 ! P=rho*Rd*T = [n*m_air*0.001 kg/g]*Rd*T
3127 ! = n*[m_air*0.001*Rd]*T
3128 ! = nRT
3129 ! Further changes required for generalization
3130 !
3131 ! 1) Conversion for surface emissions not included as yet (if any is
3132 ! needed)
3133 !
3134 implicit none
3135
3136 ! Arguments:
3137 real(8), intent(inout) :: model_col(obsoper%nmodlev) ! Model-space profile to have its units changed
3138 logical, optional, intent(in) :: ppb_opt
3139 logical, optional, intent(in) :: incr_opt
3140
3141 ! Locals:
3142 real(8) :: zcoef
3143 integer :: exp_P,exp_T
3144 logical :: ppb_out, incr_out
3145 real(8), parameter :: rho_stp=1.293 ! kg/m^3
3146
3147 if (obsoper%constituentId < 0) return
3148
3149 ! No conversion necessary for these BUFR numbers
3150 if (any( obsoper%varno == (/ BUFR_UNIT_OptDepth,BUFR_UNIT_OptDepth2, &
3151 BUFR_UNIT_OptDepth3, BUFR_UNIT_MR_NVaerosol, BUFR_NETT /) )) return
3152
3153 if (present(ppb_opt)) then
3154 ppb_out = ppb_opt
3155 else
3156 ppb_out = .false.
3157 end if
3158
3159 if (present(incr_opt)) then
3160 incr_out = incr_opt
3161 else
3162 incr_out = .false.
3163 end if
3164
3165 zcoef = 1.
3166 exp_T = 0 ! exponent of multiplicative factor obsoper%tt
3167 exp_P = 0 ! exponent of multiplicative factor obsoper%pp
3168
3169 ! Convert to ug/kg if not already in those units
3170
3171 if (obsoper%varName(1:2) == 'AF' .or. obsoper%varName(1:2) == 'AC') then
3172
3173 ! PM2.5 or PM10
3174
3175 if (any(obsoper%varno == (/ BUFR_UNIT_VMR, BUFR_UNIT_VMR2, &
3176 BUFR_UNIT_MolePerMole, BUFR_UNIT_MolePerMole2, &
3177 BUFR_UNIT_NumberDensity, BUFR_UNIT_MolarDensity, &
3178 BUFR_UNIT_PartPress, BUFR_UNIT_PartPress2, &
3179 BUFR_UNIT_DU, BUFR_UNIT_DU2, BUFR_UNIT_DU3, BUFR_UNIT_DU4, &
3180 BUFR_UNIT_IntegND, BUFR_UNIT_IntegND2 /) )) then
3181
3182 call utl_abort("oopc_convertUnits: BUFR # " // trim(utl_str(obsoper%varno)) // " is not valid for PM" )
3183 end if
3184
3185 ! Conversion from ug/m^3 to ug/kg (scaling by Rd*T/P)
3186 zcoef = zcoef * MPC_RGAS_DRY_AIR_R8
3187 exp_T = exp_T+1 ! multiply by T
3188 exp_P = exp_P-1 ! divide by P
3189
3190 else if (obsoper%varName(1:2) == 'HU') then
3191
3192 if (.not.incr_out) then
3193 ! Converts specific humidity (q) to mass mixing ratio mmr = q/(1-q)
3194 model_col = model_col / (1.0d0 - model_col)
3195 else
3196 ! For conversion of q increment (dq) to mass mixing ratio increment (dmmr)
3197 ! dmmr = dq/(1-q)^2
3198 model_col = model_col/(1.0d0 - obsoper%trial)**2
3199 end if
3200 ! Conversion factor for kg/kg to ug/kg
3201 zcoef = zcoef * 1.0d9
3202
3203 end if
3204
3205 ! Convert from ug/kg to desired unit if ppb_out = .false.
3206
3207 if (.not.ppb_out) then
3208 select case (obsoper%varno)
3209
3210 ! The first four cases below are for integral observations which
3211 ! require a conversion, in this routine, from ug/kg to the units of
3212 ! the integrand values for integrals in pressure (Pascal). Comment C above.
3213 ! Note: 1 ug/kg = 1 ppb = 1E9 mmr (mass mixing ratio)
3214
3215 case(BUFR_UNIT_IntegDens, BUFR_UNIT_IntegDens2, BUFR_UNIT_IntegDens3)
3216
3217 ! For conversion from ug/kg to integrand values in kg/(m^2*Pa)
3218 ! Note: 1 kg/(m^2**Pa) = = 1 mmr / RG = 1E-9 ug/kg / RG
3219 !
3220 zcoef = zcoef * 1.0d-9 / ec_rg
3221
3222 case(BUFR_UNIT_IntegMolarDens)
3223
3224 ! For conversion from ug/kg to integrand values in moles/(m^2*Pa)
3225 ! Note: 1 moles/(m^2*Pa) = 1E-9 ug/kg / [RG * (1E-3 kg/g * m_gas)]
3226 !
3227 ! To convert from kg/m^2 for the gas to moles/m^2, one must
3228 ! divide by the molar mass of the gas in kg/mole.
3229 !
3230 ! Note: One u or Da (unified atomic mass unit or dalton) is numerically equivalent to 1 g/mole.
3231 ! So 1 kg is equivalent to (1E3/(atomic mass)) moles
3232
3233 zcoef = zcoef * 1.0d-6 / (ec_rg*vnl_varMassFromVarNum(obsoper%constituentId))
3234
3235 case(BUFR_UNIT_IntegND, BUFR_UNIT_IntegND2)
3236
3237 ! For conversion from ug/kg to integrand values in molecules/(m^2*Pa)
3238 ! Note: 1 molecule/(m^2*Pa) = 1E-9 ug/kg * Na / [RG * (1E-3 kg/g * m_gas)]
3239 !
3240 ! To convert from kg/m^2 for the gas to molecules/m^2, one must
3241 ! divide by the gas molar mass (kg/mole) and multiply by the Avogrado number
3242
3243 zcoef = zcoef * 1.0d-6 * MPC_AVOGADRO_R8 &
3244 / (ec_rg*vnl_varMassFromVarNum(obsoper%constituentId))
3245 case(BUFR_UNIT_DU, BUFR_UNIT_DU2, BUFR_UNIT_DU3, BUFR_UNIT_DU4)
3246
3247 ! For conversion from ug/kg to integrand values in DU/Pa
3248 ! Note: 1 DU/Pa = 1E-9 ug/kg / [RG * m_gas * rho_stp/m_air * 1E-5 m]
3249 !
3250 ! 1 DU = 0.01 mm of gas at STP = 1E-5 m of gas at STP
3251 ! = integral of gas number density at STP over 1E-5 m.
3252 ! = Na*P/(RT) * 1E-5 at STP (Na=Avogadro's number)
3253 ! = Na*(molar density) * 1E-5 at STP
3254 ! = Na*rho(STP)/(molar mass) * 1E-5 m
3255 ! = integral of air number density at STP over 1E-5 m.
3256 ! = Na*rho(air,STP)/m_air * 1.E-5 m
3257 !
3258 ! Hence 1 DU equivalent to Na*rho_stp/m_air * 1E-5 m (= 2.69E20 molecules/m^2)
3259 !
3260 ! To convert from kg/m^2 for the gas to molecules/m^2, one must
3261 ! divide by the gas molar mass (kg/mole) and multiply by the Avogrado number
3262 !
3263 ! To convert from molecules/m^2 to DU, one must divide by 2.69E20 or (Na*rho_stp/m_air * 1E-5).
3264 ! So for conversion from kg/m^2 to DU, must divide by (m_gas*rho_stp/m_air * 1E-5)
3265
3266 zcoef = zcoef * 1.0d-4 * MPC_MOLAR_MASS_DRY_AIR_R8 &
3267 /(vnl_varMassFromVarNum(obsoper%constituentId)*ec_rg*rho_stp)
3268
3269 case(BUFR_UNIT_Density, BUFR_UNIT_Density2, &
3270 BUFR_UNIT_AirDensity, BUFR_UNIT_PMDensity)
3271
3272 ! For conversion from ug/kg to kg/m^3
3273 !
3274 ! rho(gas) = mass mixing ratio * rho(air) = mass mixing ratio * P/Rd/T
3275
3276 zcoef = zcoef * 1.0d-9 / MPC_RGAS_DRY_AIR_R8
3277 exp_T = exp_T-1 ! divide by T
3278 exp_P = exp_P+1 ! multiply by P
3279
3280 case(BUFR_UNIT_MMR, BUFR_UNIT_MMR2)
3281
3282 ! For conversion from ug/kg to kg/kg
3283
3284 zcoef = zcoef * 1.0d-9
3285
3286 case(BUFR_UNIT_PartPress, BUFR_UNIT_PartPress2)
3287
3288 ! For conversion from ug/kg to partial pressure (PA)
3289 !
3290 ! parial pressure = P * vmr
3291 ! = P * m_air/m_gas * mass mixing ratio
3292
3293 zcoef = zcoef * 1.0d-9 * MPC_MOLAR_MASS_DRY_AIR_R8 &
3294 /vnl_varMassFromVarNum(obsoper%constituentId)
3295 exp_P = exp_P+1 ! multiply by P
3296
3297 case(BUFR_UNIT_NumberDensity)
3298
3299 ! For conversion from ug/kg to molecules/m^3
3300 !
3301 ! Number density of gas = Na*rho(gas)/m_gas = Na*rho(air) * mass mixing ratio /m_gas
3302 ! = Na * P/Rd/T * mass mixing ratio /m_gas
3303
3304 zcoef = zcoef * 1.0d-6 * MPC_AVOGADRO_R8/MPC_RGAS_DRY_AIR_R8 &
3305 /vnl_varMassFromVarNum(obsoper%constituentId)
3306 exp_T = exp_T-1 ! divide by T
3307 exp_P = exp_P+1 ! multiply by P
3308
3309 case(BUFR_UNIT_MolarDensity)
3310
3311 ! For conversion from ug/kg to moles/m^3
3312 !
3313 ! Mole density of gas = rho(gas)/m_gas = rho(air) * mass mixing ratio /m_gas
3314 ! = P/Rd/T * mass mixing ratio /m_gas
3315
3316 zcoef = zcoef * 1.0d-6 /MPC_RGAS_DRY_AIR_R8 &
3317 /vnl_varMassFromVarNum(obsoper%constituentId)
3318 exp_T = exp_T-1 ! divide by T
3319 exp_P = exp_P+1 ! multiply by P
3320
3321 case(BUFR_UNIT_VMR, BUFR_UNIT_VMR2, BUFR_UNIT_MolePerMole, &
3322 BUFR_UNIT_MolePerMole2)
3323
3324 ! For conversion from ug/kg to vmr (or moles/mole)
3325
3326 zcoef = zcoef * 1.0d-9 * MPC_MOLAR_MASS_DRY_AIR_R8 &
3327 /vnl_varMassFromVarNum(obsoper%constituentId)
3328
3329 case default
3330
3331 call utl_abort('oopc_convertUnits: Unknown obs units ' // &
3332 'for varno = ' // trim(utl_str(obsoper%varno)) )
3333
3334 end select
3335 end if
3336
3337 ! Apply constant scaling
3338 model_col = model_col * zcoef
3339
3340 if (exp_T /= 0) then
3341 if (any(obsoper%tt <= 0.)) then
3342 call utl_abort("oopc_convertUnits: " // &
3343 "Missing valid temperature for conversion.")
3344 end if
3345 model_col = model_col * obsoper%tt**exp_T
3346 end if
3347
3348 if (exp_P /= 0) model_col = model_col * obsoper%pp**exp_P
3349
3350 end subroutine oopc_convertUnits
3351
3352 !--------------------------------------------------------------------------
3353 ! oopc_required_field
3354 !--------------------------------------------------------------------------
3355 function oopc_required_field(varName,varno) result(needed)
3356 !
3357 !:Purpose: To determine whether the specifed field name is required
3358 ! somewhere in the observation operators for a particular
3359 ! observation type.
3360 !
3361 implicit none
3362
3363 ! Arguments:
3364 character(len=*), intent(in) :: varName ! Name of field
3365 integer, intent(in) :: varno ! BUFR descriptor element
3366 ! Result:
3367 logical :: needed
3368
3369 select case(trim(varName))
3370 case('TT')
3371
3372 select case (varno)
3373 case(BUFR_UNIT_Density,BUFR_UNIT_Density2,BUFR_UNIT_AirDensity, &
3374 BUFR_UNIT_PMDensity,BUFR_UNIT_NumberDensity,BUFR_UNIT_MolarDensity)
3375 needed = .true.
3376 case default
3377 needed = .false.
3378 end select
3379
3380 case default
3381 needed = .false.
3382 end select
3383
3384 end function oopc_required_field
3385
3386 !--------------------------------------------------------------------------
3387 ! oopc_vertObsLayersWgts
3388 !--------------------------------------------------------------------------
3389 subroutine oopc_vertObsLayersWgts(operator,ixtr,success,kmode,skipType)
3390 !
3391 !:Purpose: To calculate integration (or averaging) weights "wgts" required for vertical
3392 ! integration (or averaging) for the full vertical range or a set of target layers.
3393 ! Given the calculated weights and a user input array vector X, the integral
3394 ! (or average) for a given layer i would be given by sum(wgts(i,:)*X(:))
3395 !
3396 !:Arguments:
3397 ! :operator: Operator type: 'integ' or 'avg'
3398 ! :ixtr: Flag indicating if obs outside model vertical range: 0 for no.
3399 ! :successs: Success of integration (or averaging)
3400 ! :kmode: Observation model stage used to allow option of tropo
3401 ! increment determination from total (or avg) column data when
3402 ! obsoper%columnBound > obsoper%vlayertop for
3403 ! nobslev=1. For kmode=3, calc only for region between
3404 ! obsoper%columnBound and surface. For kmode=2, split
3405 ! calc for model top to obsoper%columnBound and
3406 ! obsoper%columnBound and surface.
3407 !
3408 ! :skipType: Skipping processing of specific target layers depending on case:
3409 ! 'default' - skipping application via input success_opt only
3410 ! 'doAll&noExtrap' - application of both success_opt and outbound_opt
3411 !
3412 implicit none
3413
3414 ! Arguments:
3415 integer, intent(inout) :: ixtr(obsoper%nobslev) ! Flag indicating if obs outside model vertical range: 0 for no.
3416 logical, intent(inout) :: success(obsoper%nobslev) ! success of integration (or averaging)
3417 integer, intent(in) :: kmode ! Mode of observation operator
3418 character(len=*), intent(in) :: operator ! Operator type: 'integ' or 'avg'
3419 character(len=*), intent(in) :: skipType ! Skipping processing of specific target layers depending on case
3420
3421 ! Locals:
3422 integer :: obslevIndex,tropo_mode
3423 real(8) :: vlayertop_ref,vlayerbottom_ref,modlevindexBot_ref,checkID
3424
3425 ! Conduct initial setup for vertical integration (or avegaging) components
3426
3427 call ppo_vertLayersSetup(operator,obsoper%pp,obsoper%nmodlev)
3428
3429 ! Ensure that each layer is within model vertical range.
3430
3431 do obslevIndex=1,obsoper%nobslev
3432
3433 if (obsoper%vlayerbottom(obslevIndex) < obsoper%vlayertop(obslevIndex)) then
3434 success(1:obsoper%nobslev)=.false.
3435 write(*,*) 'oopc_vertObsLayersWgts: WARNING. ' // &
3436 'Layer top/bot value problem.', &
3437 obsoper%vlayertop(obslevIndex), &
3438 obsoper%vlayerbottom(obslevIndex), &
3439 '. Entire profile skipped over.'
3440 return
3441 else if (obsoper%vlayerbottom(obslevIndex) < obsoper%pp(1)*1.01 .or. &
3442 obsoper%vlayertop(obslevIndex) > obsoper%pp(obsoper%nmodlev)*0.99) then
3443
3444 success(obslevIndex)=.false.
3445 if (obsoper%vlayerbottom(obslevIndex) < obsoper%pp(1)*1.01) then
3446 ixtr(obslevIndex)=1
3447 else
3448 ixtr(obslevIndex)=2
3449 end if
3450 write(*,*) 'oopc_vertObsLayersWgts: WARNING. Layer top/bot value problem.', &
3451 obsoper%vlayertop(obslevIndex), obsoper%vlayerbottom(obslevIndex)
3452 cycle
3453 end if
3454 if (obsoper%vlayerbottom(obslevIndex) > &
3455 obsoper%pp(obsoper%nmodlev)*0.999) then
3456
3457 obsoper%vlayerbottom(obslevIndex)=obsoper%pp(obsoper%nmodlev)*0.999
3458
3459 end if
3460 if (obsoper%vlayertop(obslevIndex) < obsoper%pp(1)*1.001) then
3461 obsoper%vlayertop(obslevIndex)=obsoper%pp(1)*1.001
3462 end if
3463
3464 end do
3465
3466 tropo_mode=0
3467
3468 if (obsoper%nobslev == 1 .and. kmode >= 2 .and. trim(operator) == 'integ' &
3469 .and. obsoper%vlayerbottom(1) > obsoper%pp(obsoper%nmodlev)*0.99 .and. &
3470 obsoper%constituentId >= 0 .and. obsoper%vco == 4) then
3471
3472 ! Check if constituent id is recognized (function will abort if not recognized)
3473 if ( obsoper%constituentId >= 0 ) then
3474 checkID = vnl_varMassFromVarnum(obsoper%constituentId)
3475 end if
3476
3477 vlayerbottom_ref=obsoper%vlayerbottom(1)
3478
3479 ! Check for special treatment if tropo_mode>=1, kmode=2,3, and nobslev=1 for
3480 ! column observations (obsoper%vco=4) that extend to the surface.
3481
3482 tropo_mode = oopc_tropo_mode(obsoper%constituentId)
3483
3484 if ( tropo_mode >= 1 .and. obsoper%columnBound > obsoper%vlayertop(1) ) then
3485
3486 if (obsoper%iavgkern /= 0) then
3487 call utl_abort("oopc_vertObsLayersWgts: Use of averaging ' // &
3488 'kernels not possible with reduced range of increment profile.")
3489 end if
3490
3491 if (kmode == 2 .and. tropo_mode == 1) then
3492
3493 ! When kmode=2, split calc in two. This is done due to difference in
3494 ! calc at the interface region when producing zh and zhp. The tangent linear
3495 ! model in the lower region for kmode=2 must be consistent with
3496 ! that associated with kmode=3.
3497
3498 ! Start with bottom region in order to use correct zhp with oopc_genOper
3499 ! when use of this operator is requested.
3500
3501 vlayertop_ref=obsoper%vlayertop(1)
3502
3503 obsoper%vlayertop(1) = obsoper%columnBound
3504
3505 call ppo_vertIntegWgts(obsoper%vlayertop,obsoper%vlayerbottom, &
3506 obsoper%nmodlev,obsoper%nobslev,obsoper%modlevindexTop, &
3507 obsoper%modlevindexBot,obsoper%zh,wgts_opt=obsoper%zhp, &
3508 skipType_opt=skipType,outbound_opt=ixtr,success_opt=success)
3509
3510 ! Apply generalized innovation operator if requested
3511 if (obsoper%applyGenOper) call oopc_genoper(kmode)
3512
3513 obsoper%applyGenOper=.false.
3514
3515 modlevindexBot_ref=obsoper%modlevindexBot(1)
3516
3517 ! Reset top and bottom values for integration (or averaging) of the remaining region.
3518 ! The second integration (or averaging) provides the change in upper level contributions to the
3519 ! total column (or average) from the assimilation of other observations.
3520
3521 obsoper%vlayertop(1)=vlayertop_ref
3522 obsoper%vlayerbottom(1)=obsoper%columnBound
3523 else
3524 ! Reset top new value. Restricts adjoint/tangent linear calcs to this reduced region.
3525 obsoper%vlayertop(1)=obsoper%columnBound
3526 end if
3527
3528 end if
3529 else
3530 vlayerbottom_ref=0.0
3531 end if
3532
3533 ! Calculate vertical integration (or averaging) components for specified obs layer(s).
3534
3535 if ( trim(operator) == 'integ' ) then
3536
3537 call ppo_vertIntegWgts(obsoper%vlayertop,obsoper%vlayerbottom, &
3538 obsoper%nmodlev,obsoper%nobslev,obsoper%modlevindexTop, &
3539 obsoper%modlevindexBot,obsoper%zh,wgts_opt=obsoper%zhp, &
3540 skipType_opt=skipType,outbound_opt=ixtr,success_opt=success, &
3541 dealloc_opt=.true.)
3542
3543 ! If tropo_mode=1, reset original vertical range
3544 ! for the tangent linear operator
3545 if (obsoper%nobslev == 1 .and. kmode == 2 .and. &
3546 obsoper%constituentId >= 0 .and. obsoper%vco == 4) then
3547
3548 if (tropo_mode == 1 .and. &
3549 obsoper%columnBound > obsoper%vlayertop(1) .and. &
3550 vlayerbottom_ref > obsoper%pp(obsoper%nmodlev)*0.99) then
3551
3552 obsoper%vlayerbottom(1)=vlayerbottom_ref
3553 obsoper%modlevindexBot(1)=modlevindexBot_ref
3554 end if
3555 end if
3556 else
3557
3558 call ppo_vertAvgWgts(obsoper%vlayertop,obsoper%vlayerbottom, &
3559 obsoper%nmodlev,obsoper%nobslev,obsoper%modlevindexTop, &
3560 obsoper%modlevindexBot,obsoper%zh,wgts_opt=obsoper%zhp, &
3561 skipType_opt=skipType,outbound_opt=ixtr,success_opt=success, &
3562 dealloc_opt=.true.)
3563
3564 end if
3565
3566 end subroutine oopc_vertObsLayersWgts
3567
3568 !--------------------------------------------------------------------------
3569 ! oopc_genOper
3570 !--------------------------------------------------------------------------
3571 subroutine oopc_genOper(kmode)
3572 !
3573 !:Purpose: Set generalized innovation operator for integral or layer avg
3574 ! obs. Relevant only for incremental fields. This version is
3575 ! intended to vertically distribute (approximately) the obs increments
3576 ! proportionally to the (a) background state or (b) the differences
3577 ! between a reference/climatological state and the background state.
3578 ! Option (b) is recommended for a spinup period when the shape of the
3579 ! background state is not physically representative. This is to force
3580 ! the analysis shape toward the reference/climatology shape in the absence
3581 ! of local profiles observations.
3582 !
3583 !:Input:
3584 ! :kmode: index indicating if the operator is to be applied
3585 ! :obsoper%zh,zhp: see routine ppo_vertIntegSetup
3586 ! :oopc_genOperConstraintType: index specifying the reference state
3587 ! :oopc_genOperhCorrlenExpnt: Exponent for horiz. correl. length weighting
3588 ! :oopc_genOperOmaStatsFactor: Additional OmAStats normalization factor
3589 ! :bgStats: structure containing the background stats data
3590 !
3591 !:Output:
3592 ! :obsoper%zh(obsoper%nmodlev): a*w: Final innovation model array
3593 ! (other than conversion constants)
3594 ! :obsoper%zhp(obsoper%nmodlev): w (see comments section)
3595 !
3596 ! Comments:
3597 !
3598 ! (1) This routine prepares an alternative innovation operator g,
3599 ! called the generalized innovation operator, to take the place of
3600 ! the innovation (TLM) operator h (row of zh). The operator g is
3601 ! specified as:
3602 !
3603 ! g = a*w
3604 !
3605 ! where the innovation weight operator 'w' can be set as,
3606 ! for the 1D case,
3607 !
3608 ! w = P[ (h'x)^T ] * B^{-1}
3609 !
3610 ! with h' is the part of h which excludes resolution dependence
3611 ! (only/mostly contains the physics part of h; zhp),
3612 !
3613 ! - x is the state profile rval
3614 ! - P is a window cutoff operator (sets small values to zero)
3615 ! - B is the original/initial total "vertical" covariance matrix
3616 ! (2D)
3617 !
3618 !
3619 ! and 'a' is a proportionality constant ensuring that the
3620 ! innovation increment remains unchanged for the 1D case in the
3621 ! absence of other obs., i.e.,
3622 !
3623 ! a^2 = (h*B*h^T)(w*B*w^T)^{-1},
3624 !
3625 ! Application of the state profile x (rval) is to make the
3626 ! increment profile be more proportional to the state profile.
3627 !
3628 ! The presence of B^{-1} is to negate the weight re-distribution
3629 ! from the later application of B in grad(Jo). Its specification
3630 ! is approximate to reduce the effect of numerical error in the
3631 ! generating the inverse of B.
3632 !
3633 ! A scaling by the 1/bgStats%hcorrlen^q profile is addded to
3634 ! partially mitigate effect of the variation in the vertical
3635 ! of the influence of neighbouring obs via background horizontal
3636 ! error correlations.
3637 !
3638 ! The specification of the q value is to roughly give (a) a more
3639 ! constant-like increment in the absence of the trial profile form
3640 ! constraint and (b), in the presence of the form constraint, a
3641 ! better match of the max increment level to the max concentration
3642 ! level.
3643 !
3644 ! As the horizontal density of obs is somewhat variable
3645 ! and not known locally, the specified q value is not optimal for each
3646 ! individual case but better than not applying any related mitigation.
3647 !
3648 ! (2) The matrix B is the total error covariance matrix (in physical space)
3649 ! with the related error correlation matrix
3650 ! *corvert* and std dev provided from *bCovarSetupChem_mod.ftn90*.
3651 !
3652 ! (3) NOTE: Cases with ensemble-based and or lam-based background
3653 ! covariances are not taken into account in this version.
3654 !
3655 implicit none
3656
3657 ! Arguments:
3658 integer, intent(in) :: kmode ! Index specifying if content to be applied (i.e. if kmode>1)
3659
3660 ! Locals:
3661 real(8), parameter :: pwin=0.01
3662 integer :: obslevIndex,modlevIndex,irmse,varIndex
3663 real(8) :: zwbw,zhbh,za,work(obsoper%nmodlev),fdeStddev(obsoper%nmodlev,2)
3664 real(8), parameter :: threshold=1.D-20
3665 real(8) :: zmin,rvalw(obsoper%nmodlev),rvalr(obsoper%nmodlev)
3666 real(8) :: rvalc(obsoper%nmodlev),rmse
3667 character(len=22) :: code
3668
3669 if (kmode <= 1) return
3670
3671 ! Retrieve from stored background error std dev [elemements (:,1-2)] at obs location [and inverses at elements (:,3-4)]
3672 fdeStddev = bcsc_retrieveBgStddev(obsoper%nmodlev,2,obsoper%obs_index)
3673 if (obsoper%constituentId == 0) then
3674 ! To mitigate combined effects of approximations in expressions below and
3675 ! the forced strong reduction of fdeStddev(:,1) in the top level(s) of the
3676 ! model. Bkgd ozone error std dev was set to near zero at the lid. The
3677 ! absence of a resettting such as below would otherwise
3678 ! lead to oscillatory numerical error effects in the first few levels.
3679 ! This change will still force a near-zero (small) increment at the lid.
3680 fdeStddev(1,2) = fdeStddev(2,2)
3681 end if
3682
3683 ! Identify variable position index in background error correlation matrices
3684
3685 varIndex=1
3686 do while (trim(bgStats%varNameList(varIndex)) /= '')
3687 if (trim(bgStats%varNameList(varIndex)) == trim(obsoper%varName)) exit
3688 varIndex=varIndex+1
3689 end do
3690
3691 if (trim(bgStats%varNameList(varIndex)) == '') then
3692 call utl_abort('oopc_genOper: Background stats not found for ' // &
3693 trim(obsoper%varName) )
3694 end if
3695
3696 ! Initialize reference mass (mixing ratio) weighting profile as profile from trial field
3697
3698 rvalr(1:obsoper%nmodlev)=obsoper%trial(1:obsoper%nmodlev)
3699
3700 ! Check on rvalr
3701
3702 if (any(abs(rvalr(:)) <= threshold)) then
3703 if (all(abs(rvalr(:)) <= threshold)) then
3704 rvalr(:)=1.0
3705 else
3706 zmin=minval(abs(rvalr(:)), &
3707 mask=abs(rvalr(:)) > threshold)
3708 where (abs(rvalr(:)) <= threshold) rvalr(:)=zmin
3709 end if
3710 end if
3711
3712 ! Loop over obs elements
3713
3714 write(code,'(I22)') obsoper%obs_index
3715
3716 do obslevIndex=1,obsoper%nobslev
3717
3718 if (.not.obsoper%success(obslevIndex)) cycle
3719
3720 ! Set reference mass (mixing ratio) weighting profile
3721
3722 rvalw(1:obsoper%nmodlev)=rvalr(1:obsoper%nmodlev)
3723 if (trim(oopc_genOperConstraintType(obsoper%constituentId)) == 'Diff') then
3724
3725 ! Set reference mass weighting profile according to the difference between
3726 ! an external reference (such as a climatology) and trial field profiles.
3727 !
3728 ! This is a mechanism to force the solution profile shape somewhat towards that
3729 ! of the external reference when the reliability of the vertical structure of the
3730 ! trial field is not high.
3731 !
3732 ! This can be used at the beginning of long assimilation periods if the initial
3733 ! trial field is not as realistic as may be desired or somewhat mitigate
3734 ! gradual biasing of vertical structures that might otherwise occur from assimilation
3735 ! of integrated quantities (when there is insufficient data from other observation types).
3736 !
3737 ! The larger the rms difference of xc (external reference) with xb (trial field profile),
3738 ! the greater is the influence of this difference in the weighting. If there is
3739 ! little difference, then either xb or xc can be directly used as the weighting
3740 ! profile; xb is used below.
3741
3742 rmse=0.0d0
3743 irmse=0
3744 rvalc(1:obsoper%nmodlev) = oopc_getProfile(oopc_bgRef,code)
3745
3746 zmin=pwin*maxval(abs(obsoper%zhp(obslevIndex,1:obsoper%nmodlev)))
3747 do modlevIndex=1,obsoper%nmodlev
3748 if (obsoper%zhp(obslevIndex,modlevIndex) > zmin) then
3749 irmse=irmse+1
3750 rmse=rmse+((rvalc(modlevIndex)-obsoper%trial(modlevIndex))* &
3751 fdeStddev(modlevIndex,2))**2
3752 end if
3753 end do
3754 if (irmse > 0) rmse=rmse*2.0d0/irmse
3755 if (rmse < 1.0d0) then
3756 rvalw(1:obsoper%nmodlev)=rvalc(1:obsoper%nmodlev)
3757 else
3758 rvalw(1:obsoper%nmodlev)=rvalc(1:obsoper%nmodlev)- &
3759 obsoper%trial(1:obsoper%nmodlev)
3760 where (abs(rvalw(1:obsoper%nmodlev)) < 0.01d0* &
3761 rvalc(1:obsoper%nmodlev))
3762 rvalw(1:obsoper%nmodlev)=sign(0.01d0*rvalc(1:obsoper%nmodlev), &
3763 rvalw(1:obsoper%nmodlev))
3764 end where
3765 end if
3766 if (all(abs(rvalw(:)) <= threshold)) then
3767 rvalw(1:obsoper%nmodlev)=rvalc(1:obsoper%nmodlev)
3768 else if (any(abs(rvalw(:)) <= threshold)) then
3769 zmin=minval(abs(rvalw(:)),mask=abs(rvalw(:)) > threshold)
3770 where (abs(rvalw(:)) <= threshold) rvalw(:)= zmin
3771 end if
3772 end if
3773
3774 ! Begin preparation of the new innovation operator w (=new zhp)
3775
3776 if (obsoper%nobslev == 1 .and. obsoper%modlevindexTop(obslevIndex) == 1 .and. &
3777 obsoper%modlevindexBot(obslevIndex) == obsoper%nmodlev) then
3778
3779 ! Treat as total column obs. Here, zhp would otherwise be approx. equal
3780 ! to 1 except for the near-end points of the model vertical domain,
3781 ! the latter due to the discretized domain. Not using zhp avoids this
3782 ! discretization issue from weakly affecting results at the boundaries.
3783
3784 work(1:obsoper%nmodlev)=1.0
3785
3786 else
3787
3788 ! Account for localized obs function (e.g. partial columns, Jacobians. For Jacobians,
3789 ! zhp must also be independent of the model layer thicknesses.)
3790
3791 ! Apply cutoff
3792 zmin=pwin*maxval(abs(obsoper%zhp(obslevIndex,1:obsoper%nmodlev)))
3793 where (abs(obsoper%zhp(obslevIndex,1:obsoper%nmodlev)) < zmin)
3794 work(1:obsoper%nmodlev)=0.0d0
3795 elsewhere
3796 work(1:obsoper%nmodlev)=obsoper%zhp(obslevIndex,1:obsoper%nmodlev)
3797 endwhere
3798 end if
3799
3800 ! Application of 1D space vertical covariance inverse B^{-1} for partially
3801 ! mitigate the weight impact of the later application of B in finalizing
3802 ! grad(Jo). Note: fdeStddev(:,2)=1.0/fdeStddev(:,1)
3803
3804 ! Also impose
3805 ! - approximate mitigation of horizontal correlation effect
3806 ! - shape forcing profile via rvalw
3807
3808 work(1:obsoper%nmodlev)=work(1:obsoper%nmodlev)*rvalw(1:obsoper%nmodlev) &
3809 *fdeStddev(1:obsoper%nmodlev,2)
3810
3811 obsoper%zhp(obslevIndex,:)=0.0D0
3812 !$OMP PARALLEL DO PRIVATE(modlevIndex)
3813 do modlevIndex=obsoper%modlevindexTop(obslevIndex), &
3814 obsoper%modlevindexBot(obslevIndex)
3815 obsoper%zhp(obslevIndex,modlevIndex) = fdeStddev(modlevIndex,2) &
3816 /bgStats%hcorrlen(modlevIndex,varIndex)**oopc_genOperHCorrlenExpnt(varIndex) &
3817 *sum(work(1:obsoper%nmodlev) &
3818 *bgStats%corverti(1:obsoper%nmodlev,modlevIndex,varIndex))
3819 end do
3820 !$OMP END PARALLEL DO
3821
3822 ! Determine proportionality factor 'a' = (h*B*h^T)(w*B*w^T)^{-1}
3823
3824 ! First determine/estimate w*B*w^T (zwbw)
3825
3826 !$OMP PARALLEL DO PRIVATE(modlevIndex)
3827 do modlevIndex=obsoper%modlevindexTop(obslevIndex), &
3828 obsoper%modlevindexBot(obslevIndex)
3829
3830 work(modlevIndex)=sum(obsoper%zhp(obslevIndex, &
3831 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) &
3832 *bgStats%corvert(modlevIndex,obsoper%modlevindexTop(obslevIndex): &
3833 obsoper%modlevindexBot(obslevIndex),varIndex) &
3834 *fdeStddev(obsoper%modlevindexTop(obslevIndex): &
3835 obsoper%modlevindexBot(obslevIndex),1))*fdeStddev(modlevIndex,1)
3836 end do
3837 !$OMP END PARALLEL DO
3838
3839 zwbw=sum(obsoper%zhp(obslevIndex,obsoper%modlevindexTop(obslevIndex): &
3840 obsoper%modlevindexBot(obslevIndex)) &
3841 *work(obsoper%modlevindexTop(obslevIndex): &
3842 obsoper%modlevindexBot(obslevIndex)))
3843
3844 ! Determine/estimate h*B*h^T (zhbh)
3845
3846 !$OMP PARALLEL DO PRIVATE(modlevIndex)
3847 do modlevIndex=obsoper%modlevindexTop(obslevIndex), &
3848 obsoper%modlevindexBot(obslevIndex)
3849
3850 work(modlevIndex)=sum(obsoper%zh(obslevIndex, &
3851 obsoper%modlevindexTop(obslevIndex):obsoper%modlevindexBot(obslevIndex)) &
3852 *bgStats%corvert(modlevIndex,obsoper%modlevindexTop(obslevIndex): &
3853 obsoper%modlevindexBot(obslevIndex),varIndex) &
3854 *fdeStddev(obsoper%modlevindexTop(obslevIndex): &
3855 obsoper%modlevindexBot(obslevIndex),1))*fdeStddev(modlevIndex,1)
3856 end do
3857 !$OMP END PARALLEL DO
3858
3859 zhbh=sum(obsoper%zh(obslevIndex,obsoper%modlevindexTop(obslevIndex): &
3860 obsoper%modlevindexBot(obslevIndex)) &
3861 *work(obsoper%modlevindexTop(obslevIndex): &
3862 obsoper%modlevindexBot(obslevIndex)))
3863
3864 ! Set proportionality factor 'a'
3865
3866 za=sqrt(zhbh/zwbw)*oopc_genOperOmAStatsFactor(varIndex)
3867
3868 ! Set final innovation operator
3869
3870 obsoper%zh(obslevIndex,1:obsoper%nmodlev)= &
3871 obsoper%zhp(obslevIndex,1:obsoper%nmodlev)*za
3872
3873 end do
3874
3875 end subroutine oopc_genOper
3876
3877 !--------------------------------------------------------------------------
3878 ! oopc_getColBoundary
3879 !--------------------------------------------------------------------------
3880 function oopc_getColBoundary(iconstituentId,nmodlev,pressmod,tt,height,hu_opt, &
3881 uu_opt,vv_opt) result(boundPress)
3882 !
3883 !:Purpose: To determine and to store the boundary (e.g. tropopause or PBL)
3884 ! pressure levels if needed by the observation operators.
3885 implicit none
3886
3887 ! Arguments:
3888 integer, intent(in) :: iconstituentId ! BUFR code element of Table 08046 identifying the constituent
3889 integer, intent(in) :: nmodlev ! number of model levels for variables other than uu and vv
3890 real(8), intent(in) :: pressmod(nmodlev)! model pressure array
3891 real(8), intent(in) :: tt(nmodlev) ! model temperature (Kelvin)
3892 real(8), intent(in) :: height(nmodlev) ! height (meters)
3893 real(8), optional, intent(in) :: hu_opt(nmodlev) ! specific humidity
3894 real(8), optional, intent(in) :: uu_opt(:) ! model zonal wind component(m/s)
3895 real(8), optional, intent(in) :: vv_opt(:) ! model meridional wind component (m/s)
3896 ! Result:
3897 real(8) :: boundPress ! pressure level of boundary to be imposed
3898
3899 ! Locals:
3900 integer :: tropo_bound
3901
3902 boundPress = -1.0d0
3903
3904 if (oopc_tropo_mode(iconstituentId) == 0) return
3905
3906 tropo_bound = oopc_tropo_bound(iconstituentId)
3907
3908 if (tropo_bound > 0) then
3909 if (.not.all(tt < 0.) .and. .not.all(height < 0.) ) then
3910
3911 select case(tropo_bound)
3912 case(1)
3913
3914 ! Get tropopause pressure level
3915
3916 boundPress = phf_get_tropopause(nmodlev,pressmod,tt,height,hu_opt=hu_opt)
3917
3918 case(2)
3919
3920 ! Get PBL pressure level
3921
3922 boundPress = phf_get_pbl(nmodlev,pressmod,tt,height,hu_opt=hu_opt, &
3923 uu_opt=uu_opt,vv_opt=vv_opt)
3924
3925 case default
3926 call utl_abort("oopc_getColBoundary: Unrecognized value for tropo_bound of " &
3927 // trim(utl_str(tropo_bound)) )
3928 end select
3929
3930 end if
3931 end if
3932
3933 ! Use tropo_column_top value if tropo_bound=0 or model derived boundary was unsuccessful
3934 if (boundPress < 0.0d0) boundPress = oopc_tropo_column_top(iconstituentId)
3935
3936 end function oopc_getColBoundary
3937
3938 !--------------------------------------------------------------------------
3939 ! oopc_addColBoundary
3940 !--------------------------------------------------------------------------
3941 subroutine oopc_addColBoundary(headerIndex,boundPress)
3942 !
3943 !:Purpose: To add column boundary data to oopc_columnBoundary which can be
3944 ! retrieved later using a header index.
3945 !
3946 implicit none
3947
3948 ! Arguments:
3949 integer, intent(in) :: headerIndex ! Header index
3950 real(8), intent(in) :: boundPress ! Pressure boundary
3951
3952 if (.not.associated(oopc_columnBoundary%data1d)) then
3953 call oss_obsdata_alloc(oopc_columnBoundary,oopc_obsdata_maxsize, dim1=1)
3954 oopc_columnBoundary%nrep = 0
3955 end if
3956
3957 ! In this case nrep will count the number of filled reps in the data arrays
3958 oopc_columnBoundary%nrep = oopc_columnBoundary%nrep+1
3959
3960 if (oopc_columnBoundary%nrep > oopc_obsdata_maxsize) then
3961 call utl_abort('oopc_addColBoundary: Reach max size of array ' // &
3962 trim(utl_str(oopc_obsdata_maxsize)) )
3963 end if
3964
3965 ! Use the header number as the unique code for this obs data
3966 write(oopc_columnBoundary%code(oopc_columnBoundary%nrep),'(I22)') headerIndex
3967
3968 oopc_columnBoundary%data1d(1,oopc_columnBoundary%nrep) = boundPress
3969
3970 end subroutine oopc_addColBoundary
3971
3972 !--------------------------------------------------------------------------
3973 ! oopc_retrieveColBoundary
3974 !--------------------------------------------------------------------------
3975 function oopc_retrieveColBoundary(headerIndex) result(boundPress)
3976 !
3977 !:Purpose: To retrieve previously saved column boundary data in
3978 ! oopc_columnBoundary from the header index.
3979 !
3980 implicit none
3981
3982 ! Arguments:
3983 integer, intent(in) :: headerIndex
3984 ! Result:
3985 real(8) :: boundPress
3986
3987 ! Locals:
3988 character(len=22) :: code
3989
3990 write(code,'(I22)') headerIndex
3991
3992 boundPress = oss_obsdata_get_element(oopc_columnBoundary,code,1)
3993
3994 end function oopc_retrieveColBoundary
3995
3996 !==========================================================================
3997 !-------- Stand-alone routines to read and extract climatology fields -----
3998
3999 ! Could be placed in a category 4 module.
4000 ! public :: oopc_readFields, oopc_addToProfileSet, oopc_getProfile
4001
4002 !--------------------------------------------------------------------------
4003 ! oopc_readFields
4004 !--------------------------------------------------------------------------
4005 subroutine oopc_readFields(climatFields,filename,variable, &
4006 maxNumFields,maxNumTypes, &
4007 fieldRequired,success,filetype_opt)
4008 !
4009 !:Purpose: To read climatrology (reference) fields as directed by input
4010 !
4011 ! Comments:
4012 ! - Fields are provided in RPN/fst files
4013 ! - Reference fields can be in a separate RPN file with name provided
4014 ! within 'filename' if filetype='TXT' or provided as 'filename' if it
4015 ! refers to an RPN standaard file.
4016 ! - Fields assumed to be of the same units as those of the
4017 ! corresponding input trial fields
4018 !
4019 implicit none
4020
4021 ! Arguments:
4022 type(struct_oopc_field), intent(out) :: climatFields(0:maxNumFields,maxNumTypes)
4023 character(len=*), intent(in) :: filename
4024 integer, intent(in) :: maxNumFields,maxNumTypes
4025 logical, intent(out) :: success
4026 character(len=*), intent(in) :: variable
4027 character(len=*), intent(in) :: fieldRequired(0:maxNumFields)
4028 character(len=*), optional, intent(in) :: filetype_opt
4029
4030 ! Locals:
4031 character(len=3) :: filetype
4032 character(len=256) :: fname
4033 character(len=4) :: varName
4034 character(len=12) :: etiket
4035 integer :: varIndex,id,nd,j,numvar,ijour,imonth,iday,itime,latIndex
4036 real(8) :: day
4037 integer :: datestamp
4038 integer, external :: newdate
4039 integer, external :: fnom, fclos
4040 integer :: ierr, nulun, ios
4041 logical :: fileExists
4042 logical :: timeInterp
4043 integer :: ni, nj, nkeys, kind
4044 real(8), allocatable :: array1(:,:,:),array2(:,:,:),lvls(:),xlat(:),xlong(:)
4045 real(8), allocatable :: pressclim(:),ozoneclim(:,:)
4046 character (len=128) :: ligne
4047
4048 ! Initialize dimensions to zero
4049
4050 climatFields(:,:)%nlon=0
4051 climatFields(:,:)%nlat=0
4052 climatFields(:,:)%nlev=1
4053
4054 if ( trim(variable) == 'CH' ) then
4055 if ( all(fieldRequired(:) == 'Trial') ) then
4056 ! Not needed
4057 success=.true.
4058 return
4059 end if
4060 end if
4061
4062 inquire(file=trim(filename),exist=fileExists)
4063 if ( .not.fileExists ) then
4064 write(*,*) '----------------------------------------------------'
4065 write(*,*) 'oopc_readFields: COULD NOT FIND file ' // trim(filename)
4066 write(*,*) '----------------------------------------------------'
4067 success = .false.
4068 return
4069 else
4070 success = .true.
4071 end if
4072
4073 ! Check for file names containing climatological fields or input directives
4074
4075 if ( present(filetype_opt) ) then
4076 filetype = trim(fileType_opt)
4077 else
4078 filetype = 'RPN'
4079 end if
4080
4081 nulun=0
4082 ierr=0
4083 if ( filetype == 'TXT' ) then
4084 ierr=fnom(nulun,trim(filename),'SEQ',0)
4085 if ( ierr == 0 ) then
4086 open(unit=nulun, file=trim(filename), status='OLD')
4087 ios=0
4088
4089 if ( trim(variable) == 'CH' ) then
4090
4091 ! CH variable kind (for constituent fields)
4092
4093 read(nulun,'(A)',iostat=ios,err=10,end=10) ligne
4094 do while (trim(adjustl(ligne(1:14))) /= 'SECTION IV:')
4095 read(nulun,'(A)',iostat=ios,err=10,end=11) ligne
4096 end do
4097
4098 ! Read number of constituents with associated input file(s)
4099
4100 read(nulun,*,iostat=ios,err=10,end=10) numvar
4101 if (numvar <= 0) go to 10
4102 else
4103 numvar=1
4104 nd=1
4105 call utl_abort('oopc_readFields: Variable kind or name ' // &
4106 trim(variable) // ' not taken into account')
4107 end if
4108 end if
4109 else if ( filetype == 'RPN' ) then
4110 numvar=1
4111 nd=1
4112 else if ( filetype /= 'RPN' ) then
4113 call utl_abort('oopc_readFields: File type ' // trim(filetype) // &
4114 ' not recognized')
4115 end if
4116
4117 if ( ierr /= 0 ) then
4118 call utl_abort('oopc_readFields: COULD NOT OPEN file ' // trim(filename))
4119 end if
4120
4121 ! Initialization
4122
4123 timeInterp = .true.
4124 datestamp=tim_getDateStamp()
4125 ierr = newdate(datestamp,ijour,itime,-3)
4126 if ( ierr < 0 ) then
4127 call utl_abort('oopc_readFields: Invalid datestamp ' // &
4128 trim(utl_str(datestamp)) )
4129 end if
4130 imonth = MOD(ijour/100,100)
4131 iday = MOD(ijour,100)
4132 day=iday+itime*1.0D-8
4133 if (day > 15.) then
4134 day=day-15.0
4135 else
4136 day=day+15.0
4137 end if
4138
4139 ! Get needed fields for each file/varIndex
4140
4141 do varIndex=1,numvar
4142
4143 if ( trim(variable) == 'CH' ) then
4144
4145 ! Read id,nd
4146 ! id: constituent identifier code; (0 for ozone, ...)
4147 ! nd: number of sets; 1 or 2 (nd=2 required when different profile
4148 ! sets need to be merged according to the tropopause height
4149 ! when the first set referring to strato files and teh second
4150 ! to tropo fields)
4151
4152 read(nulun,*,iostat=ios,err=10,end=10)
4153 read(nulun,*,iostat=ios,err=10,end=10) id,nd
4154 varName=vnl_varnameFromVarnum(0,id)
4155
4156 read(nulun,*,iostat=ios,err=10,end=10) fname
4157 inquire(file=trim(fname),exist=fileExists)
4158 if ( .not. fileExists ) then
4159 call utl_abort('oopc_readFields: Did not find file ' // trim(fname))
4160 end if
4161 else
4162 id=varIndex
4163 ! Currently assumes nunmar = 1 and fname = filename. Could be extended
4164 fname = filename
4165 varname = trim(variable)
4166 end if
4167
4168 do j=1,nd
4169
4170 if ( trim(fname) == 'ozoneclim98' ) then
4171 timeInterp = .false.
4172 call ozo_read_Climatology(datestamp,nlat_opt=nj,nlev_opt=nkeys, &
4173 press_opt=pressclim,ozone_opt=ozoneclim)
4174 id=0
4175 ni=1
4176 allocate(array1(1,nj,nkeys),lvls(nkeys),xlat(nj),xlong(1))
4177 ! Convert from ppmv to microgram/kg
4178 array1(1,1:nj,1:nkeys) = ozoneclim(1:nj,1:nkeys) * &
4179 MPC_MOLAR_MASS_O3_R8 / (1.0d-3 * MPC_MOLAR_MASS_DRY_AIR_R8)
4180 lvls(1:nkeys) = pressclim(1:nkeys)
4181 deallocate(ozoneclim,pressclim)
4182 kind = 2
4183 xlong(1)=0.0d0
4184 do latIndex = 1, nj
4185 xlat(latIndex) = (latIndex-1)*180.0d0/(nj-1) - 90.0d0
4186 end do
4187 etiket = ' '
4188 else
4189 if ( nd == 2 ) then
4190 read(nulun,*,iostat=ios,err=10,end=10) etiket
4191 else
4192 etiket = ' '
4193 end if
4194 call utl_readFstField(trim(fname),varName,-1,imonth,-1,etiket, &
4195 ni,nj,nkeys,array1,xlat_opt=xlat,xlong_opt=xlong, &
4196 lvls_opt=lvls,kind_opt=kind)
4197 end if
4198
4199 climatFields(id,j)%nlon=ni
4200 climatFields(id,j)%nlat=nj
4201 climatFields(id,j)%nlev=nkeys
4202 climatFields(id,j)%ivkind=kind
4203
4204 allocate(climatFields(id,j)%field(ni,nj,nkeys))
4205 allocate(climatFields(id,j)%vlev(nkeys),climatFields(id,j)%lon(ni))
4206 allocate(climatFields(id,j)%lat(nj))
4207
4208 climatFields(id,j)%lat(1:nj)=xlat(1:nj)*MPC_RADIANS_PER_DEGREE_R8
4209 climatFields(id,j)%lon(1:ni)=xlong(1:ni)*MPC_RADIANS_PER_DEGREE_R8
4210 where (climatFields(id,j)%lon(1:ni) < 0.0)
4211 climatFields(id,j)%lon(1:ni)=2.0*MPC_PI_R8 + climatFields(id,j)%lon(1:ni)
4212 end where
4213 climatFields(id,j)%vlev(1:nkeys)=lvls(1:nkeys)
4214
4215 if (.not.timeInterp) then
4216
4217 climatFields(id,j)%field(:,:,:) = array1(:,:,:)
4218
4219 else
4220
4221 ! Following for interpolation as a function of days from mid-months.
4222
4223 if (iday > 15) then
4224 if (imonth == 12) then
4225 call utl_readFstField(trim(fname),varName,-1,1,-1,etiket, &
4226 ni,nj,nkeys,array2,lvls_opt=lvls,kind_opt=kind)
4227 else
4228 call utl_readFstField(trim(fname),varName,-1,imonth+1,-1, &
4229 etiket,ni,nj,nkeys,array2,lvls_opt=lvls,kind_opt=kind)
4230 end if
4231
4232 ! Linearly interpolate in time
4233 ! (approximately - assumes 30 day months)
4234
4235 climatFields(id,j)%field(:,:,:) = (array1(:,:,:)*(30.0-day)+array2(:,:,:)*day)/30.0
4236
4237 else if (iday <= 15) then
4238 if (imonth == 1) then
4239 call utl_readFstField(trim(fname),varName,-1,12,-1,etiket, &
4240 ni,nj,nkeys,array2,lvls_opt=lvls,kind_opt=kind)
4241 else
4242 call utl_readFstField(trim(fname),varName,-1,imonth-1,-1, &
4243 etiket,ni,nj,nkeys,array2,lvls_opt=lvls,kind_opt=kind)
4244 end if
4245
4246 ! Linearly interpolate in time
4247 ! (approximately - assumes 30 day months)
4248
4249 climatFields(id,j)%field(:,:,:) = (array2(:,:,:)* &
4250 (30.0-day)+array1(:,:,:)*day)/30.0
4251 end if
4252 end if
4253
4254 if (allocated(array1)) deallocate(array1,lvls,xlat,xlong)
4255 if (allocated(array2)) deallocate(array2)
4256
4257 end do
4258 end do
4259
4260 10 if (ios > 0) then
4261 call utl_abort('oopc_readFields: READING PROBLEM.' // &
4262 ' File read error message number: ' // trim(utl_str(ios)))
4263 end if
4264 close(unit=nulun)
4265 ierr = fclos(nulun)
4266 if ( any(fieldRequired(:) == 'Diff') .and. trim(variable) == 'CH' ) then
4267 do j=0,maxNumFields
4268 if ( climatFields(j,1)%nlon == 0 .and. trim(fieldRequired(j)) == 'Diff' ) then
4269 call utl_abort('oopc_readFields: READING PROBLEM. Did not' // &
4270 ' find SECTION IV required for constituent ID ' // &
4271 trim(utl_str(j)))
4272 end if
4273 end do
4274 end if
4275
4276 return
4277
4278 11 close(unit=nulun)
4279 ierr = fclos(nulun)
4280 if ( any(fieldRequired(:) == 'Diff') .and. trim(variable) == 'CH' ) then
4281 call utl_abort('oopc_readFields: READING PROBLEM. Did not find ' // &
4282 'SECTION IV.')
4283 end if
4284
4285 end subroutine oopc_readFields
4286
4287 !--------------------------------------------------------------------------
4288 ! oopc_addToProfileSet
4289 !--------------------------------------------------------------------------
4290 subroutine oopc_addToProfileSet(climatFields,climatProfileSet,maxNumFields,maxNumTypes, &
4291 numModelLevs,modelPressLevs,modelHeightLevs,obsLat, &
4292 obsLong,obsIndex,maxsize,varKind_opt,varNumber_opt,tt_opt,hu_opt)
4293
4294 !:Purpose: To determine and to store a profile at obs location as part of a cumulative
4295 ! profile set for a specific variable
4296 !
4297 !:Input:
4298 !
4299 ! :climatFields: Input fields from which interpolations are done
4300 ! :climatProfileSet: Input profile set
4301 ! :maxNumFields: Size of first dimension for climatFields
4302 ! :maxNumTypes: Size of second dimension for climatFields
4303 ! :numModelLevs: Number of model levels
4304 ! :modelPressLevs Model pressure array (Pa)
4305 ! :modelHeightLevs: Model height (m)
4306 ! :obsLat: Latitude (rad)
4307 ! :obsLong: Longitude (rad)
4308 ! :obsIndex: Unique measurement identifier
4309 ! :varKind_opt: variable kind (currently only relevant for 'CH')
4310 ! :varNumber_opt: Constituent id
4311 ! :tt_opt: Model temperature (Kelvin)
4312 ! :hu_opt: Specific humidity
4313 ! :maxsize: Max number of obs for which climatProfileSet will be used
4314 !
4315 !:Output:
4316 !
4317 ! :climatProfileSet: Updated profile set (with one profile added for (obs_long,obs_lat))
4318 !
4319 implicit none
4320
4321 ! Arguments:
4322 type(struct_oopc_field), intent(in) :: climatFields(0:maxNumFields,maxNumTypes)
4323 type(struct_oss_obsdata), intent(inout) :: climatProfileSet
4324 integer, intent(in) :: maxNumFields
4325 integer, intent(in) :: maxNumTypes
4326 integer, intent(in) :: obsIndex
4327 integer, intent(in) :: numModelLevs
4328 integer, intent(in) :: maxsize
4329 real(8), intent(in) :: modelPressLevs(numModelLevs)
4330 real(8), intent(in) :: modelHeightLevs(numModelLevs)
4331 real(8), intent(in) :: obsLat
4332 real(8), intent(in) :: obsLong
4333 integer, optional, intent(in) :: varNumber_opt
4334 real(8), optional, intent(in) :: tt_opt(:)
4335 real(8), optional, intent(in) :: hu_opt(:)
4336 character(len=*), optional, intent(in) :: varKind_opt
4337
4338 ! Locals
4339 integer :: level,start,id
4340 real(8) :: tropo_press, refprof(numModelLevs),refprof2(numModelLevs),dt
4341 real(8), allocatable :: pressrefin(:)
4342 logical, allocatable :: success(:)
4343
4344 if ( present(varNumber_opt) ) then
4345 if ( varNumber_opt < 0 ) return
4346 id = varNumber_opt
4347 else
4348 id = 0
4349 end if
4350
4351 if ( present(varKind_opt) ) then
4352 ! Not currently used
4353 end if
4354
4355 if (climatFields(id,1)%nlat == 0) return
4356
4357 ! Set vertical levels of reference.
4358 ! Convert to pressure coordinate if needed.
4359
4360 if (allocated(pressrefin)) deallocate(pressrefin)
4361 allocate(pressrefin(climatFields(id,1)%nlev))
4362 pressrefin(:) = climatFields(id,1)%vlev(1:climatFields(id,1)%nlev)
4363
4364 if (allocated(success)) deallocate(success)
4365 allocate(success(climatFields(id,1)%nlev))
4366 success(:)=.true.
4367
4368 if (climatFields(id,1)%ivkind == 2) then
4369 pressrefin(:)=pressrefin(:)*100. ! Conversion from hPa to Pa.
4370 else if (climatFields(id,1)%ivkind == 0) then
4371 where (pressrefin < modelHeightLevs(numModelLevs))
4372 pressrefin=modelHeightLevs(numModelLevs)
4373 end where
4374 pressrefin(:) = phf_convert_z_to_pressure(pressrefin,modelHeightLevs,modelPressLevs, &
4375 climatFields(id,1)%nlev,numModelLevs,obsLat,success)
4376 else if (climatFields(id,1)%ivkind == 4) then
4377 pressrefin(:)=pressrefin(:) + modelHeightLevs(numModelLevs)
4378 pressrefin(:) = phf_convert_z_to_pressure(pressrefin,modelHeightLevs,modelPressLevs, &
4379 climatFields(id,1)%nlev,numModelLevs,obsLat,success)
4380 else if (climatFields(id,1)%ivkind == 1) then
4381 pressrefin(:)=pressrefin(:)*modelPressLevs(numModelLevs) ! Convert from sigma to Pa
4382 else
4383 call utl_abort('oopc_addToProfileSet: Cannot handle vertical coordinate of kind ' // trim(utl_str(climatFields(id,1)%ivkind)))
4384 end if
4385
4386 ! Interpolate to obs lat/long (or lat) location and model level
4387
4388 call oopc_column_hbilin(climatFields(id,1)%field,pressrefin, &
4389 climatFields(id,1)%nlon,climatFields(id,1)%nlat,climatFields(id,1)%nlev, &
4390 climatFields(id,1)%lon,climatFields(id,1)%lat,obsLong,obsLat, &
4391 refprof,modelPressLevs,numModelLevs)
4392
4393 if (climatFields(id,2)%nlat > 0 .and. climatFields(id,2)%nlon > 0 &
4394 .and. climatFields(id,2)%nlev > 0) then
4395
4396 if ( .not. present(tt_opt) ) then
4397 call utl_abort('oopc_addToProfileSet: Missing TT for determining ' // &
4398 'tropopause pressure')
4399 end if
4400 if ( any(tt_opt <= 0.0d0) ) then
4401 call utl_abort('oopc_addToProfileSet: Invalid TT for determining ' // &
4402 'tropopause pressure')
4403 end if
4404
4405 ! Get second reference field (for troposphere)
4406
4407 tropo_press=-1.0
4408
4409 if ( present(hu_opt) ) then
4410 if (all(hu_opt >= 0.0D0)) then
4411 tropo_press=phf_get_tropopause(numModelLevs,modelPressLevs, &
4412 tt_opt,modelHeightLevs,hu_opt=hu_opt)
4413 else
4414 tropo_press=phf_get_tropopause(numModelLevs,modelPressLevs, &
4415 tt_opt,modelHeightLevs)
4416 end if
4417 else
4418 tropo_press=phf_get_tropopause(numModelLevs,modelPressLevs,tt_opt,modelHeightLevs)
4419 end if
4420
4421 if (tropo_press > 0) then
4422
4423 ! Set vertical levels of reference.
4424 ! Convert to pressure coordinate if needed
4425
4426 if (allocated(pressrefin)) deallocate(pressrefin)
4427 allocate(pressrefin(climatFields(id,2)%nlev))
4428 pressrefin(:)= climatFields(id,2)%vlev(1:climatFields(id,2)%nlev)
4429
4430 if (allocated(success)) deallocate(success)
4431 allocate(success(climatFields(id,2)%nlev))
4432 success(:)=.true.
4433
4434 if (climatFields(id,2)%ivkind == 2) then
4435 pressrefin(:)=pressrefin(:)*100. ! Conversion from hPa to Pa.
4436 else if (climatFields(id,2)%ivkind == 0) then
4437 where (pressrefin < modelHeightLevs(numModelLevs))
4438 pressrefin=modelHeightLevs(numModelLevs)
4439 end where
4440 pressrefin(:) = phf_convert_z_to_pressure(pressrefin, &
4441 modelHeightLevs,modelPressLevs, &
4442 climatFields(id,2)%nlev,numModelLevs, &
4443 obsLat,success)
4444 else if (climatFields(id,2)%ivkind == 4) then
4445 pressrefin(:)=pressrefin(:) + modelHeightLevs(numModelLevs)
4446 pressrefin(:) = phf_convert_z_to_pressure(pressrefin, &
4447 modelHeightLevs,modelPressLevs, &
4448 climatFields(id,2)%nlev,numModelLevs,obsLat, &
4449 success)
4450 else if (climatFields(id,2)%ivkind == 1) then
4451 pressrefin(:)=pressrefin(:)*modelPressLevs(numModelLevs) ! Convert from sigma to Pa
4452 else
4453 call utl_abort('oopc_addToProfileSet: Cannot handle vertical ' // &
4454 'coordinate of kind ' // trim(utl_str(climatFields(id,2)%ivkind)))
4455 end if
4456
4457 ! Interpolate to obs lat/long (or lat) and model levels
4458
4459 call oopc_column_hbilin(climatFields(id,2)%field,pressrefin, &
4460 climatFields(id,2)%nlon,climatFields(id,2)%nlat,climatFields(id,2)%nlev, &
4461 climatFields(id,2)%lon,climatFields(id,2)%lat,obsLong,obsLat, &
4462 refprof2,modelPressLevs,numModelLevs)
4463
4464 end if
4465
4466 ! Combine with upper level profile
4467
4468 do level=numModelLevs,3,-1
4469 if (modelPressLevs(level) < tropo_press) exit
4470 refprof(level)=refprof2(level)
4471 end do
4472 start=level
4473
4474 ! Apply linear combination of four levels just above the tropopause
4475
4476 do level=start,max(2,start-3),-1
4477 dt=(start+1.0-level)/5.0
4478 refprof(level)=dt*refprof2(level) + (1.0-dt)*refprof(level)
4479 end do
4480
4481 end if
4482
4483 if (allocated(pressrefin)) deallocate(pressrefin)
4484 if (allocated(success)) deallocate(success)
4485
4486 ! ------- Save in climatProfileSet ---------
4487
4488 if (.not.associated(climatProfileSet%data1d)) then
4489 call oss_obsdata_alloc(climatProfileSet, maxsize, dim1=numModelLevs)
4490 climatProfileSet%nrep = 0
4491 end if
4492
4493 ! Here, nrep will count the number of filled elements in the data arrays
4494 climatProfileSet%nrep = climatProfileSet%nrep+1
4495
4496 if (climatProfileSet%nrep > maxsize) then
4497 call utl_abort('oopc_addToProfilesSet: Reach max size of array ' // &
4498 trim(utl_str(maxsize)) )
4499 end if
4500
4501 ! obsIndex serves as the unique locator code
4502 write(climatProfileSet%code(climatProfileSet%nrep),'(I22)') obsIndex
4503
4504 ! Save profile in climatProfileSet
4505
4506 climatProfileSet%data1d(:,climatProfileSet%nrep) = refprof(:)
4507
4508 end subroutine oopc_addToProfileSet
4509
4510 !--------------------------------------------------------------------------
4511 ! oopc_getProfile
4512 !--------------------------------------------------------------------------
4513 function oopc_getProfile(climatProfileSet,code) result(profile)
4514 !
4515 !:Purpose: To extract and provide profile from climatProfileSet according to
4516 ! code value.
4517 !
4518 implicit none
4519
4520 ! Arguments
4521 type(struct_oss_obsdata), intent(inout) :: climatProfileSet ! Profile set
4522 character(len=*), intent(in) :: code ! unique obs identifying code
4523 ! Result:
4524 real(8) :: profile(climatProfileSet%dim1) ! retrieved array from obsdata%data1d of dimension obsdata%dim1
4525
4526 ! Locals:
4527 integer :: status ! search success (0 = found; 1 = no data; 2 = not found)
4528
4529 profile = oss_obsdata_get_array1d(climatProfileSet,code,status)
4530 if (status > 0) then
4531 call utl_abort("oopc_getProfile: Code not found - " // trim(code))
4532 end if
4533
4534 end function oopc_getProfile
4535
4536 !--------------------------------------------------------------------------
4537 ! oopc_column_hbilin
4538 !--------------------------------------------------------------------------
4539 subroutine oopc_column_hbilin(field,vlev,nlong,nlat,nlev,xlong,xlat, &
4540 plong,plat,vprof,vlevout,nlevout)
4541 !
4542 !:Purpose: Horizontal bilinear interpolation from a 3D field to a profile at (plong,plat).
4543 ! Assumes vertical interpolation not needed or already done.
4544 !
4545 ! This version can be used with fields that are not part of the background state,
4546 ! such as climatologies.
4547 !
4548 ! This version does not depend in column_data and gridstatevector modules.
4549 !
4550 implicit none
4551
4552 ! Arguments:
4553 integer, intent(in) :: nlong ! number or longitudes
4554 integer, intent(in) :: nlat ! number or latitudes
4555 integer, intent(in) :: nlev ! number of vertical levels
4556 integer, intent(in) :: nlevout ! number of target vertical levels
4557 real(8), intent(in) :: field(nlong,nlat,nlev) ! 3D field
4558 real(8), intent(in) :: vlev(nlev) ! vertical levels of input field (in pressure)
4559 real(8), intent(in) :: xlong(nlong) ! longitudes (radians)
4560 real(8), intent(in) :: xlat(nlat) ! latitudes (radians)
4561 real(8), intent(in) :: plong ! target longitude (radians)
4562 real(8), intent(in) :: plat ! target latitude (radian)
4563 real(8), intent(in) :: vlevout(nlevout) ! target vertical levels (in pressure)
4564 real(8), intent(out) :: vprof(nlevout) ! profile at (plong,plat)
4565
4566 ! Locals:
4567 real(8) :: lnvlev(nlev),lnvlevout(nlevout),plong2
4568 integer :: ilev,lonIndex,latIndex,i,j
4569 real(8) :: DLDX, DLDY, DLDP, DLW1, DLW2, DLW3, DLW4
4570
4571 call utl_tmg_start(30,'--StateToColumn')
4572
4573 ! Find near lat/long grid points
4574
4575 if ( nlong > 1 ) then
4576 plong2 = plong
4577 if (plong2 < 0.0) plong2 = 2.D0*MPC_PI_R8 + plong2
4578 do lonIndex = 2, nlong
4579 if (xlong(lonIndex-1) < xlong(lonIndex)) then
4580 if (plong2 >= xlong(lonIndex-1) .and. plong2 <= xlong(lonIndex)) exit
4581 else
4582 ! Assumes this is a transition between 360 to 0 (if it exists). Skip over.
4583 end if
4584 end do
4585 lonIndex = lonIndex-1
4586 else
4587 lonIndex=0
4588 end if
4589
4590 do latIndex = 2, nlat
4591 if (plat <= xlat(latIndex)) exit
4592 end do
4593 latIndex = latIndex-1
4594
4595 if ( lonIndex == 0 ) then
4596
4597 ! Set lat interpolation weights
4598
4599 DLDY = (plat - xlat(latIndex))/(xlat(latIndex+1)-xlat(latIndex))
4600
4601 DLW1 = (1.d0-DLDY)
4602 DLW2 = DLDY
4603
4604 ! Set vertical interpolation weights (assumes pressure vertical coordinate)
4605
4606 lnvlevout(:) = log(vlevout(:))
4607 lnvlev(:) = log(vlev(:))
4608
4609 ilev = 1
4610 do i = 1, nlevout
4611 do j = ilev, nlev
4612 if (lnvlevout(i) < lnvlev(j)) exit ! assumes lnvlevout and lnvlev increase with index
4613 end do
4614 ilev = j-1
4615 if (ilev < 1) then
4616 ilev = 1
4617 else if (ilev >= nlev) then
4618 ilev = nlev-1
4619 end if
4620
4621 DLDP = (lnvlev(ilev+1)-lnvlevout(i))/(lnvlev(ilev+1)-lnvlev(ilev))
4622
4623 vprof(i) = DLDP* (DLW1 * field(lonIndex,latIndex,ilev) &
4624 + DLW2 * field(lonIndex,latIndex+1,ilev)) &
4625 + (1.d0-DLDP)* (DLW1 * field(lonIndex,latIndex,ilev+1) &
4626 + DLW2 * field(lonIndex,latIndex+1,ilev+1))
4627 end do
4628
4629 else
4630
4631 ! Set lat/long interpolation weights
4632
4633 DLDX = (plong - xlong(lonIndex))/(xlong(lonIndex+1)-xlong(lonIndex))
4634 DLDY = (plat - xlat(latIndex))/(xlat(latIndex+1)-xlat(latIndex))
4635
4636 DLW1 = (1.d0-DLDX) * (1.d0-DLDY)
4637 DLW2 = DLDX * (1.d0-DLDY)
4638 DLW3 = (1.d0-DLDX) * DLDY
4639 DLW4 = DLDX * DLDY
4640
4641 ! Set vertical interpolation weights (assumes pressure vertical coordinate)
4642
4643 lnvlevout(:) = log(vlevout(:))
4644 lnvlev(:) = log(vlev(:))
4645
4646 ilev = 1
4647 do i = 1, nlevout
4648 do j = ilev, nlev
4649 if (lnvlevout(i) < lnvlev(j)) exit ! assumes lnvlevout and lnvlev increase with index
4650 end do
4651 ilev = j-1
4652 if (ilev < 1) then
4653 ilev = 1
4654 else if (ilev >= nlev) then
4655 ilev = nlev-1
4656 end if
4657
4658 DLDP = (lnvlev(ilev+1)-lnvlevout(i))/(lnvlev(ilev+1)-lnvlev(ilev))
4659
4660 vprof(i) = DLDP* (DLW1 * field(lonIndex,latIndex,ilev) &
4661 + DLW2 * field(lonIndex+1,latIndex,ilev) &
4662 + DLW3 * field(lonIndex,latIndex+1,ilev) &
4663 + DLW4 * field(lonIndex+1,latIndex+1,ilev)) &
4664 + (1.d0-DLDP)* (DLW1 * field(lonIndex,latIndex,ilev+1) &
4665 + DLW2 * field(lonIndex+1,latIndex,ilev+1) &
4666 + DLW3 * field(lonIndex,latIndex+1,ilev+1) &
4667 + DLW4 * field(lonIndex+1,latIndex+1,ilev+1))
4668 end do
4669 end if
4670
4671 call utl_tmg_stop(30)
4672
4673 end subroutine oopc_column_hbilin
4674
4675end module obsOperatorsChem_mod
