-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsetuprad.f90
2671 lines (2357 loc) · 121 KB
/
setuprad.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module rad_setup
implicit none
private
public:: setup
interface setup; module procedure setuprad; end interface
contains
subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,&
obstype,isis,is,rad_diagsave,init_pass,last_pass)
!$$$ subprogram documentation block
! . . . .
! subprogram: setuprad compute rhs of oi equation for radiances
! prgmmr: derber org: np23 date: 1995-07-06
!
! abstract: read in data, first guess, and obtain rhs of oi equation
! for radiances.
!
! program history log:
! 1995-07-06 derber
! 1996-11-xx wu, data from prepbufr file
! 1996-12-xx mcnally, changes for diagnostic file and bugfix
! 1998-04-30 weiyu yang mpi version
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2003-12-23 kleist - remove sigma assumptions (use pressure)
! 2004-05-28 kleist - subroutine call update
! 2004-06-17 treadon - update documenation
! 2004-07-23 weng,yan,okamoto - incorporate MW land and snow/ice emissivity
! models for AMSU-A/B and SSM/I
! 2004-08-02 treadon - add only to module use, add intent in/out
! 2004-10-06 parrish - modifications for nonlinear qc
! 2004-10-15 derber - modify parts of IR quality control
! 2004-10-28 treadon - replace parameter tiny with tiny_r_kind
! 2004-11-22 derber - remove weight, add logical for boundary point
! 2004-11-30 xu li - add SST physical retrieval algorithm
! 2004-12-22 treadon - add outer loop number to name of diagnostic file
! 2005-01-20 okamoto - add ssm/i radiance assimilation
! 2005-01-22 okamoto - add TB jacobian with respect to ocean surface wind
! through MW ocean emissivity model
! 2005-02-22 derber - alter surface determination and improve quality control
! 2005-02-28 treadon - increase size of character variable holding diagnostic
! file name
! 2005-03-02 derber - modify use of surface flages and quality control
! and adjoint of surface emissivity
! 2005-03-04 xu li - restructure code related to sst retrieval
! 2005-03-07 todling,treadon - place lower bound on sum2
! 2005-03-09 parrish - nonlinear qc change to account for inflated obs error
! 2005-03-16 derber - save observation time
! 2005-04-11 treadon - add logical to toggle on/off nonlinear qc code
! 2005-04-18 treadon - modify sections of code related to sst retrieval
! 2005-06-01 treadon - add code to load/use extended vertical profile arrays in rtm
! 2005-07-06 derber - modify for mhs and hirs/4
! 2005-07-29 treadon - modify tnoise initialization; add varinv_use
! 2005-09-20 xu,pawlak - modify sections of code related to ssmis
! 2005-09-28 derber - modify for new radinfo and surface info input from read routines
! 2005-10-14 derber - input grid location and fix regional lat/lon
! 2005-10-17 treadon - generalize accessing of elements from obs array
! 2005-10-20 kazumori - modify sections of code related to amsre
! 2005-11-04 derber - place lower bound (0.0) on computed clw
! 2005-11-14 li - modify avhrr related code
! 2005-11-18 treadon - correct thin snow test to apply to microwave
! 2005-11-18 kazumori - modify sections of amsre diagnostic file
! 2005-11-29 parrish - remove call to deter_sfc_reg (earlier patch for regional mode)
! 2005-12-16 derber - add check on skin temperature to clw bias correction
! 2005-12-20 derber - add transmittance qc check to mw sensors
! 2006-01-09 treadon - introduce get_ij
! 2006-01-12 treadon - replace pCRTM with CRTM
! 2006-01-31 todling - add obs time to output diag files
! 2006-02-01 liu - add ssu
! 2006-02-02 treadon - rename prsi(l) as ges_prsi(l)
! 2006-02-03 derber - add new obs control and change printed stats
! 2006-03-21 treadon - add optional perturbation to observation
! 2006-03-24 treadon - bug fix - add iuse_rad to microwave channel varinv check
! 2006-04-19 treadon - rename emisjac as dtbduv_on (accessible via obsmod)
! 2006-04-27 derber - remove rad_tran_k, process data one profile at a time
! write data in jppf chunks
! 2006-05-10 derber - add check on maximum number of levels for RT
! 2006-05-30 derber,treadon - modify diagnostic output
! 2006-06-06 su - move to wgtlim to constants module
! 2006-07-27 kazumori - modify factor of bc predictor(clw) for AMSR-E
! and input of qcssmi
! 2006-07-28 derber - modify to use new inner loop obs data structure
! - unify NL qc and add satellite and solar azimuth angles
! 2006-07-31 kleist - change call to intrppx, no longer get ps at ob location
! 2006-12-21 sienkiewicz - add 'no85GHz' flag for F8 SSM/I
! 2007-01-24 kazumori- modify to qcssmi subroutine output and use ret_ssmis
! for ssmis_las only (assumed UKMO SSMIS data)
! 2007-03-09 su - remove the perturbation to the observation
! 2007-03-19 tremolet - binning of observations
! 2007-04-04 wu - do not load ozone jacobian if running regional mode
! 2007-05-30 h.liu - replace c1 with constoz in ozone jacobian
! 2007-06-05 tremolet - add observation diagnostics structure
! 2007-06-08 kleist/treadon - add prefix (task id or path) to diag_rad_file
! 2007-06-29 jung - update CRTM interface
! 2008-01-30 h.liu/treadon - add SSU cell pressure correction block
! 2008-05-21 safford - rm unused vars and uses
! 2008-12-03 todling - changed handle of tail%time
! 2009-12-07 b.yan - changed qc for channel 5 (relaxed)
! 2009-08-19 guo - changed for multi-pass setup with dtime_check(), and
! new arguments init_pass and last_pass.
! 2009-12-08 guo - cleaned diag output rewind with open(position='rewind')
! - fixed a bug in diag header output while is not init_pass.
! 2010-03-01 gayno - allow assimilation of "mixed" amsua fovs
! 2010-03-30 collard - changes for interface with CRTM v2.0.
! 2010-03-30 collard - Add CO2 interface (fixed value for now).
! 2010-04-08 h.liu -add SEVIRI assimilation
! 2010-04-16 hou/kistler add interface to module ncepgfs_ghg
! 2010-04-29 zhu - add option newpc4pred for new preconditioning for predictors
! 2010-05-06 zhu - add option adp_anglebc variational angle bias correction
! 2010-05-13 zhu - add option passive_bc for bias correction of passive channels
! 2010-05-19 todling - revisit intrppx CO2 handle
! 2010-06-10 todling - reduce pointer check by getting CO2 pointer at this level
! - start adding hooks of aerosols influence on RTM
! 2010-07-15 kleist - reintroduce capability to write out predictor terms (not predicted bias) and
! pressure level that corresponds to peak of weighting function
! 2010-07-16 yan - update quality control of mw water vapor sounding channels (amsu-b and mhs)
! - add a new input (tbc) to in call qcssmi(..) and
! remove 'ssmis_uas,ssmis_las,ssmis_env,ssmis_img' in call qcssmi(..)
! Purpose: to keep the consistent changes with qcssmi.f90
! 2010-08-10 wu - setup corresponding vegetation types (nmm_to_crtm) for IGBP in regional
! parameter nvege_type: old=24, IGBP=20
! 2010-08-17 derber - move setup input and crtm call to crtm_interface (intrppx) to simplify routine
! 2010-09-30 zhu - re-order predterms and predbias
! 2010-12-16 treadon - move cbias update before calc_clw
! 2011-02-17 todling - add knob to turn off O3 Jacobian from IR instruments (per Emily Liu's work)
! 2011-03-13 li - (1) associate nst_gsi and nstinfo (use radinfo) to handle nst fields
! - (2) modify to save nst analysis related diagnostic variables
! 2011-04-07 todling - newpc4pred now in radinfo
! 2011-05-04 todling - partially merge in Min-Jeong Kim's cloud clear assimilation changes (connect to Metguess)
! 2011-05-16 todling - generalize handling of jacobian matrix entries
! 2011-05-20 mccarty - updated for ATMS
! 2011-06-08 zhu - move assignments of tnoise_cld values to satinfo file via varch_cld, use lcw4crtm
! 2011-06-09 sienkiewicz - call to qc_ssu needs tb_obs instead of tbc
! 2011-07-10 zhu - add jacobian assignments for regional cloudy radiance
! 2011-09-28 collard - Fix error trapping for CRTM failures.
! 2012-05-12 todling - revisit opts in gsi_metguess_get (4crtm)
! 2012-11-02 collard - Use cloud detection channel flag for IR.
! 2013-02-13 eliu - Add options for SSMIS instruments
! - Add two additional bias predictors for SSMIS radiances
! - Tighten up QC checks for SSMIS
! 2013-02-19 sienkiewicz - add adjustable preweighting for SSMIS bias terms
! 2013-07-10 zhu - add upd_pred as an update indicator for bias correction coeficitient
! 2013-07-19 zhu - add emissivity sensitivity predictor for radiance bias correction
! 2013-11-19 sienkiewicz - merge back in changes for adjustable preweighting for SSMIS bias terms
! 2013-11-21 todling - inquire diag-file version using get_radiag
! 2013-12-10 zhu - apply bias correction to tb_obs for ret_amsua calculation
! 2013-12-21 eliu - add amsu-a obs errors for allsky condition
! 2013-12-21 eliu - add error handling for CLWP calculation for allsky
! 2014-01-17 zhu - add cld_rbc_idx for bias correction sample to handle cases with cloud
! inconsistency between obs and first guess for all-sky microwave radiance
! 2014-01-19 zhu - add scattering index calculation, add it as a predictor for allsky
! - calculate retrieved clw using bias-corrected tsim
! 2014-01-28 todling - write sensitivity slot indicator (ioff) to header of diagfile
! 2014-01-31 mkim - Remove abs(60.0degree) boundary which existed for all-sky MW radiance DA
! 2014-02-01 mkim - Move all-sky mw obserr to subroutine obserr_allsky_mw
! 2014-02-05 todling - Remove overload of diagbufr slot (not allowed)
! 2014-04-17 todling - Implement inter-channel ob correlated covariance capability
! 2014-04-27 eliu - change qc_amsua/atms interface
! 2014-04-27 eliu - change call_crtm interface to output clear-sky Tb under all-sky condition (optional)
! 2014-04-27 eliu - add cloud effect calculation for AMSU-A/ATMS under all-sky condition
! 2014-05-29 thomas - add lsingleradob capability (originally of mccarty)
! 2014-08-01 zhu - remove scattering index predictor
! - add all-sky obs error adjustment based on scattering index, diff of clw,
! cloud mismatch info, and surface wind speed
! 2014-12-30 derber - Modify for possibility of not using obsdiag
! 2015-01-15 zhu - change amsua quality control interface to apply emissivity sensitivity
! screen to all-sky AMSUA and ATMS radiance
! 2015-01-16 ejones - Added call to qc_gmi for gmi observations
! - Added saphir
! 2015-02-12 ejones - Write gwp to diag file for GMI
! 2015-03-11 ejones - Added call to qc_amsr2 for amsr2 observations
! 2015-03-23 ejones - Added call to qc_saphir for saphir observations
! 2015-03-23 zaizhong ma - add Himawari-8 ahi
! 2014-08-06 todling - Correlated obs now platform-instrument specific
! 2014-09-02 todling - Must protect NST-related diag out for when NST is on
! 2014-09-03 j.jin - Added GMI 1CR radiance, obstype=gmi.
! 2014-12-30 derber - Modify for possibility of not using obsdiag
! 2015-03-31 zhu - move cloudy AMSUA radiance observation error adjustment to qcmod.f90;
! change quality control interface for AMSUA and ATMS.
! 2015-04-01 W. Gu - add isis to obs type
! 2015-08-18 W. Gu - include the dependence of the correlated obs errors on the surface types.
! 2015-09-04 J.Jung - Added mods for CrIS full spectral resolution (FSR).
! 2015-09-10 zhu - generalize enabling all-sky and aerosol usage in radiance assimilation.
! Use radiance_obstype_search & type extentions from radiance_mod.
! - special obs error & bias correction handlings are called from centralized module
! 2015-09-30 ejones - Pull AMSR2 sun azimuth and sun zenith angles for passing to quality control,
! modify qc_amsr2 function call
! 2015-10-01 guo - full res obvsr: index to allow redistribution of obsdiags
! 2016-02-15 zhu - remove the code forcing zero Jacobians for qr,qs,qg,qh for regional, let users decide
! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting
! 2016-06-24 guo - fixed the default value of obsdiags(:,:)%tail%luse to luse(n)
! . removed (%dlat,%dlon) debris.
! 2016-12-09 mccarty - add netcdf_diag capability
! 2016-07-19 W. Gu - add isis to obs type
! 2016-07-19 W. Gu - include the dependence of the correlated obs errors on the surface types
! 2016-07-19 kbathmann -move eigendecomposition for correlated obs here
! 2016-11-29 shlyaeva - save linearized H(x) for EnKF
! 2016-12-09 mccarty - add netcdf_diag capability
! 2016-03-02 mkim - Added all-sky GMI microwave radiance DA related interfaces
! 2016-03-02 mkim - Moved NCEP's AMSU-A all-sky obs error codes to a separate subroutine 'obserr_allsky_mw'
! to allow different obserror models(or coefficients) for different sensors
! while keeping setuprad.f90 compact.
! 2016-03-02 mkim - unified naming for retrieved cloud liquid water path as 'clw_obs' and obsolete
! 'clwp_amsua' and 'clw' to consider all-sky DA for other microwave sensors
! while keeping setuprad.f90 compact
! 2016-10-23 zhu - add cloudy radiance assimilation for ATMS
! 2017-02-09 guo - Remove m_alloc, n_alloc.
! . Remove my_node with corrected typecast().
! 2017-07-27 kbathmann/W. Gu -introduce rinvdiag into the rstats computation for correlated error
! 2018-04-04 zhu - add additional radiance_ex_obserr and radiance_ex_biascor calls for all-sky
! 2018-08-08 mkim - merging NCEP all-sky generalization stuff (radiance_mod, cloudy_radiance_info.txt...)
! 2018-07-24 W. Gu - Store the R-covariance matrix only needed for method=1 or 2
! 2018-07-27 W. Gu - code changes to reduce the round-off errors
! 2019-03-13 eliu - add components to handle precipitation-affected radiances
! 2019-03-13 eliu - add calculation of scattering index for MHS/ATMS
! 2019-03-27 h. liu - add ABI assimilation
!
! input argument list:
! lunin - unit from which to read radiance (brightness temperature, tb) obs
! mype - mpi task id
! nchanl - number of channels per obs
! nreal - number of pieces of non-tb information per obs
! nobs - number of tb observations to process
! obstype - type of tb observation
! isis - sensor/instrument/satellite id ex.amsua_n15
! is - integer counter for number of observation types to process
! rad_diagsave - logical to switch on diagnostic output (.false.=no output)
! channelinfo - structure containing satellite sensor information
!
! output argument list:
! aivals - array holding sums for various statistics as a function of obs type
! stats - array holding sums for various statistics as a function of channel
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use mpeu_util, only: die,perr,getindex
use kinds, only: r_kind,r_single,i_kind
use crtm_spccoeff, only: sc
use radinfo, only: nuchan,tlapmean,predx,cbias,ermax_rad,tzr_qc,&
npred,jpch_rad,varch,varch_cld,iuse_rad,icld_det,nusis,fbias,retrieval,b_rad,pg_rad,&
air_rad,ang_rad,adp_anglebc,angord,ssmis_precond,emiss_bc,upd_pred, &
passive_bc,ostats,rstats,newpc4pred,radjacnames,radjacindxs,nsigradjac,nvarjac, &
varch_sea,varch_land,varch_ice,varch_snow,varch_mixed
use gsi_nstcouplermod, only: nstinfo
use read_diag, only: get_radiag,ireal_radiag,ipchan_radiag
use guess_grids, only: sfcmod_gfs,sfcmod_mm5,comp_fact10
use m_prad, only: radheadm
use m_obsdiagNode, only: obs_diag
use m_obsdiagNode, only: obs_diags
use m_obsdiagNode, only: fptr_obsdiagNode
use m_obsdiagNode, only: obsdiagLList_nextNode
use m_obsdiagNode, only: obsdiagNode_set
use m_obsdiagNode, only: obsdiagNode_get
use m_obsdiagNode, only: obsdiagNode_assert
use obsmod, only: ianldate,ndat,mype_diaghdr,nchan_total, &
dplat,dtbduv_on,lobsdiag_forenkf,&
lobsdiagsave,nobskeep,lobsdiag_allocated,&
dirname,time_offset,lwrite_predterms,lwrite_peakwt,reduce_diag
use m_obsNode, only: obsNode
use m_radNode, only: radNode
use m_radNode, only: radNode_appendto
use m_obsLList, only: obsLList
use obsmod, only: luse_obsdiag,dval_use
use obsmod, only: netcdf_diag, binary_diag, dirname
use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, &
nc_diag_write, nc_diag_data2d, nc_diag_chaninfo_dim_set, nc_diag_chaninfo
use gsi_4dvar, only: nobs_bins,hr_obsbin,l4dvar
use gridmod, only: nsig,regional,get_ij
use satthin, only: super_val1
use constants, only: quarter,half,tiny_r_kind,zero,one,deg2rad,rad2deg,one_tenth, &
two,three,cg_term,wgtlim,r100,r10,r0_01,r_missing
use jfunc, only: jiter,miter,jiterstart
use sst_retrieval, only: setup_sst_retrieval,avhrr_sst_retrieval,&
finish_sst_retrieval,spline_cub
use m_dtime, only: dtime_setup, dtime_check
use crtm_interface, only: init_crtm,call_crtm,destroy_crtm,sensorindex,surface, &
itime,ilon,ilat,ilzen_ang,ilazi_ang,iscan_ang,iscan_pos,iszen_ang,isazi_ang, &
ifrac_sea,ifrac_lnd,ifrac_ice,ifrac_sno,itsavg, &
izz,idomsfc,isfcr,iff10,ilone,ilate, &
isst_hires,isst_navy,idata_type,iclr_sky,itref,idtw,idtc,itz_tr
use crtm_interface, only: ilzen_ang2,iscan_ang2,iszen_ang2,isazi_ang2
use clw_mod, only: calc_clw, ret_amsua, gmi_37pol_diff
use qcmod, only: qc_ssmi,qc_seviri,qc_abi,qc_ssu,qc_avhrr,qc_goesimg,qc_msu,qc_irsnd,qc_amsua,qc_mhs,qc_atms
use qcmod, only: igood_qc,ifail_gross_qc,ifail_interchan_qc,ifail_crtm_qc,ifail_satinfo_qc,qc_noirjaco3,ifail_cloud_qc
use qcmod, only: ifail_cao_qc,cao_check
use qcmod, only: ifail_iland_det, ifail_isnow_det, ifail_iice_det, ifail_iwater_det, ifail_imix_det, &
ifail_iomg_det, ifail_isst_det, ifail_itopo_det,ifail_iwndspeed_det
use qcmod, only: qc_gmi,qc_saphir,qc_amsr2
use radinfo, only: iland_det, isnow_det, iwater_det, imix_det, iice_det, &
iomg_det, itopo_det, isst_det,iwndspeed_det
use qcmod, only: setup_tzr_qc,ifail_scanedge_qc,ifail_outside_range
use state_vectors, only: svars3d, levels, svars2d, ns3d, nsdim
use oneobmod, only: lsingleradob,obchan,oblat,oblon,oneob_type
use correlated_obsmod, only: corr_adjust_jacobian, idnames
use radiance_mod, only: rad_obs_type,radiance_obstype_search,radiance_ex_obserr,radiance_ex_biascor
use sparsearr, only: sparr2, new, writearray, size, fullarray
use radiance_mod, only: radiance_ex_obserr_gmi,radiance_ex_biascor_gmi
implicit none
! Declare passed variables
type(obsLList ),target,dimension(:),intent(in):: obsLL
type(obs_diags),target,dimension(:),intent(in):: odiagLL
logical ,intent(in ) :: rad_diagsave
character(10) ,intent(in ) :: obstype
character(20) ,intent(in ) :: isis
integer(i_kind) ,intent(in ) :: lunin,mype,nchanl,nreal,nobs,is
real(r_kind),dimension(40,ndat) ,intent(inout) :: aivals
real(r_kind),dimension(7,jpch_rad),intent(inout) :: stats
logical ,intent(in ) :: init_pass,last_pass ! state of "setup" processing
! Declare external calls for code analysis
external:: stop2
! Declare local parameters
real(r_kind),parameter:: r1e10=1.0e10_r_kind
character(len=*),parameter:: myname="setuprad"
! Declare local variables
character(128) diag_rad_file
integer(i_kind) iextra,jextra,error_status
integer(i_kind) ich9,isli,icc,iccm,mm1,ixx
integer(i_kind) m,mm,jc,j,k,i
integer(i_kind) n,nlev,kval,ibin,ioff,ioff0,iii,ijacob
integer(i_kind) ii,jj,idiag,inewpc,nchanl_diag
integer(i_kind) nadir,kraintype,ierrret
integer(i_kind) ioz,ius,ivs,iwrmype
integer(i_kind) iversion_radiag, istatus
integer(i_kind) cor_opt,iinstr,chan_count
character(len=80) covtype
real(r_single) freq4,pol4,wave4,varch4,tlap4
real(r_kind) node
real(r_kind) term,tlap,tb_obsbc1,tb_obsbc16,tb_obsbc17
real(r_kind) drad,dradnob,varrad,error,errinv,useflag
real(r_kind) cg_rad,wgross,wnotgross,wgt,arg,exp_arg
real(r_kind) tzbgr,tsavg5,trop5,pangs,cld,cldp
real(r_kind) cenlon,cenlat,slats,slons,zsges,zasat,dtime
! real(r_kind) wltm1,wltm2,wltm3
real(r_kind) ys_bias_sst,cosza,val_obs
real(r_kind) sstnv,sstcu,sstph,dtp_avh,dta,dqa
real(r_kind) bearaz,sun_zenith,sun_azimuth
! real(r_kind) sfc_speed,frac_sea,clw,tpwc,sgagl,clwp_amsua,tpwc_guess_retrieval
! real(r_kind) sfc_speed,frac_sea,clw,tpwc,sgagl,tpwc_guess_retrieval
real(r_kind) sfc_speed,frac_sea,tpwc_obs,sgagl,tpwc_guess_retrieval
real(r_kind) gwp,clw_obs
real(r_kind) scat,scatp
real(r_kind) dtsavg,r90,coscon,sincon
real(r_kind) bias
real(r_kind) factch6
real(r_kind) stability,tcwv,hwp_ratio
real(r_kind) si_obs,si_fg,si_mean
logical cao_flag
logical hirs2,msu,goessndr,hirs3,hirs4,hirs,amsua,amsub,airs,hsb,goes_img,ahi,mhs,abi
type(sparr2) :: dhx_dx
real(r_single), dimension(nsdim) :: dhx_dx_array
logical avhrr,avhrr_navy,lextra,ssu,iasi,cris,seviri,atms
logical ssmi,ssmis,amsre,amsre_low,amsre_mid,amsre_hig,amsr2,gmi,saphir
logical ssmis_las,ssmis_uas,ssmis_env,ssmis_img
logical sea,mixed,land,ice,snow,toss,l_may_be_passive,eff_area
logical microwave, microwave_low
logical no85GHz
logical in_curbin, in_anybin, save_jacobian
logical account_for_corr_obs
logical,dimension(nobs):: zero_irjaco3_pole
! Declare local arrays
real(r_single),dimension(ireal_radiag):: diagbuf
real(r_single),allocatable,dimension(:,:):: diagbufex
real(r_single),allocatable,dimension(:,:):: diagbufchan
real(r_kind),dimension(npred+2):: predterms
real(r_kind),dimension(npred+2,nchanl):: predbias
real(r_kind),dimension(npred,nchanl):: pred,predchan
real(r_kind),dimension(nchanl):: err2,tbc0,raterr2,wgtjo
real(r_kind),dimension(nchanl):: varinv0
real(r_kind),dimension(nchanl):: varinv,varinv_use,error0,errf,errf0
real(r_kind),dimension(nchanl):: tb_obs,tbc,tbcnob,tlapchn,tb_obs_sdv
real(r_kind),dimension(nchanl):: tnoise,tnoise_cld
real(r_kind),dimension(nchanl):: emissivity,ts,emissivity_k
real(r_kind),dimension(nchanl):: tsim,wavenumber,tsim_bc
! real(r_kind),dimension(nchanl):: tsim_clr,tsim_clr_bc,cldeff_obs,cldeff_sim
real(r_kind),dimension(nchanl):: tsim_clr,tsim_clr_bc,cldeff_obs,cldeff_fg
real(r_kind),dimension(nsig,nchanl):: wmix,temp,ptau5
real(r_kind),dimension(nsigradjac,nchanl):: jacobian
real(r_kind),dimension(nreal+nchanl,nobs)::data_s
real(r_kind),dimension(nsig):: qvp,tvp
real(r_kind),dimension(nsig):: prsltmp
real(r_kind),dimension(nsig+1):: prsitmp
real(r_kind),dimension(nchanl):: weightmax
real(r_kind),dimension(nchanl):: cld_rbc_idx,cld_rbc_idx2
real(r_kind),dimension(nchanl):: tcc
real(r_kind) :: ptau5deriv, ptau5derivmax
real(r_kind) :: clw_guess,clw_guess_retrieval,ciw_guess,rain_guess,snow_guess,clw_avg
real(r_kind) :: tnoise_save
real(r_kind),dimension(:), allocatable :: rsqrtinv
real(r_kind),dimension(:), allocatable :: rinvdiag
!for GMI (dual scan angles)
real(r_kind),dimension(nchanl):: emissivity2,ts2, emissivity_k2,tsim2
real(r_kind),dimension(nchanl):: tsim_clr2
real(r_kind),dimension(5) :: gmi_low_angles
real(r_kind),dimension(nsig,nchanl):: wmix2,temp2,ptau52
real(r_kind),dimension(nsigradjac,nchanl):: jacobian2
real(r_kind) cosza2
integer(i_kind),dimension(nchanl):: ich,id_qc,ich_diag
integer(i_kind),dimension(nchanl):: kmax
integer(i_kind),allocatable,dimension(:) :: sc_index
integer(i_kind) :: state_ind, nind, nnz
logical channel_passive
logical,dimension(nobs):: luse
integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID
integer(i_kind):: nperobs
character(10) filex
character(12) string
type(radNode),pointer:: my_head,my_headm
type(obs_diag),pointer:: my_diag
type(obs_diags),pointer:: my_diagLL
type(rad_obs_type) :: radmod
type(obsLList),pointer,dimension(:):: radhead
type(fptr_obsdiagNode),dimension(nchanl):: odiags
logical:: muse_ii
! Notations in use: for a single obs. or a single obs. type
! nchanl : a known channel count of a given type obs stream
! nchanl_diag : a subset of "iuse"
! icc, iii : a subset of "(varinv(i)>tiny_r_kind) .and. iuse)" or qc-passed
! And for all instruments
! jpch_rad : sum(nchanl)
! nchanl_total : subset of jpch_rad, sum(icc)
radhead => obsLL(:)
save_jacobian = rad_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf
if (save_jacobian) then
ijacob = 1 ! flag to indicate jacobian saved in diagnostic file
else
ijacob = 0
endif
!**************************************************************************************
! Initialize variables and constants.
mm1 = mype+1
r90 = 90._r_kind
coscon = cos( (r90-55.0_r_kind)*deg2rad )
sincon = sin( (r90-55.0_r_kind)*deg2rad )
factch6 = zero
cld = zero
cldp = zero
tpwc_obs = zero
sgagl = zero
dtp_avh=zero
icc = 0
iccm = 0
ich9 = min(9,nchanl)
do i=1,nchanl
do j=1,npred
pred(j,i)=zero
end do
end do
! Initialize logical flags for satellite platform
cao_flag = .false.
hirs2 = obstype == 'hirs2'
hirs3 = obstype == 'hirs3'
hirs4 = obstype == 'hirs4'
hirs = hirs2 .or. hirs3 .or. hirs4
msu = obstype == 'msu'
ssu = obstype == 'ssu'
goessndr = obstype == 'sndr' .or. obstype == 'sndrd1' .or. &
obstype == 'sndrd2'.or. obstype == 'sndrd3' .or. &
obstype == 'sndrd4'
amsua = obstype == 'amsua'
amsub = obstype == 'amsub'
mhs = obstype == 'mhs'
airs = obstype == 'airs'
hsb = obstype == 'hsb'
goes_img = obstype == 'goes_img'
ahi = obstype == 'ahi'
avhrr = obstype == 'avhrr'
avhrr_navy = obstype == 'avhrr_navy'
ssmi = obstype == 'ssmi'
amsre_low = obstype == 'amsre_low'
amsre_mid = obstype == 'amsre_mid'
amsre_hig = obstype == 'amsre_hig'
amsre = amsre_low .or. amsre_mid .or. amsre_hig
amsr2 = obstype == 'amsr2'
gmi = obstype == 'gmi'
ssmis = obstype == 'ssmis'
ssmis_las = obstype == 'ssmis_las'
ssmis_uas = obstype == 'ssmis_uas'
ssmis_img = obstype == 'ssmis_img'
ssmis_env = obstype == 'ssmis_env'
iasi = obstype == 'iasi'
cris = obstype == 'cris' .or. obstype == 'cris-fsr'
seviri = obstype == 'seviri'
atms = obstype == 'atms'
saphir = obstype == 'saphir'
abi = obstype == 'abi'
ssmis=ssmis_las.or.ssmis_uas.or.ssmis_img.or.ssmis_env.or.ssmis
microwave=amsua .or. amsub .or. mhs .or. msu .or. hsb .or. &
ssmi .or. ssmis .or. amsre .or. atms .or. &
amsr2 .or. gmi .or. saphir
microwave_low =amsua .or. msu .or. ssmi .or. ssmis .or. amsre
! Determine cloud & aerosol usages in radiance assimilation
call radiance_obstype_search(obstype,radmod)
! Initialize channel related information
tnoise = r1e10
tnoise_cld = r1e10
l_may_be_passive = .false.
toss = .true.
jc=0
do j=1,jpch_rad
if(isis == nusis(j))then
jc=jc+1
if(jc > nchanl)then
write(6,*)'SETUPRAD: ***ERROR*** in channel numbers, jc,nchanl=',jc,nchanl,&
' ***STOP IN SETUPRAD***'
call stop2(71)
end if
! Load channel numbers into local array based on satellite type
ich(jc)=j
do i=1,npred
predchan(i,jc)=predx(i,j)
end do
!
! Set error instrument channels
tnoise(jc)=varch(j)
channel_passive=iuse_rad(j)==-1 .or. iuse_rad(j)==0
if (iuse_rad(j)< -1 .or. (channel_passive .and. &
.not.rad_diagsave)) tnoise(jc)=r1e10
if (passive_bc .and. channel_passive) tnoise(jc)=varch(j)
if (iuse_rad(j)>0) l_may_be_passive=.true.
if (tnoise(jc) < 1.e4_r_kind) toss = .false.
tnoise_cld(jc)=varch_cld(j)
if (iuse_rad(j)< -1 .or. (iuse_rad(j) == -1 .and. &
.not.rad_diagsave)) tnoise_cld(jc)=r1e10
if (passive_bc .and. (iuse_rad(j)==-1)) tnoise_cld(jc)=varch_cld(j)
end if
end do
if(nchanl > jc) write(6,*)'SETUPRAD: channel number reduced for ', &
obstype,nchanl,' --> ',jc
if(jc == 0 .or. toss)then
if(jc == 0 .and. mype == 0) then
write(6,*)'SETUPRAD: No channels found for ', obstype,isis
end if
if (toss .and. mype == 0) then
write(6,*)'SETUPRAD: all obs var > 1e4. do not use ',&
'data from satellite is=',isis
endif
if(nobs >0)read(lunin)
return
endif
if ( mype == 0 .and. .not.l_may_be_passive) write(6,*)mype,'setuprad: passive obs',is,isis
! Logic to turn off print of reading coefficients if not first interation or not mype_diaghdr or not init_pass
iwrmype=-99
if(mype==mype_diaghdr(is) .and. init_pass .and. jiterstart == jiter)iwrmype = mype_diaghdr(is)
! Initialize radiative transfer and pointers to values in data_s
call init_crtm(init_pass,iwrmype,mype,nchanl,nreal,isis,obstype,radmod)
! Get indexes of variables in jacobian to handle exceptions down below
ioz =getindex(radjacnames,'oz')
if(ioz>0) then
ioz=radjacindxs(ioz)
endif
ius =getindex(radjacnames,'u')
ivs =getindex(radjacnames,'v')
if(ius>0.and.ivs>0) then
ius=radjacindxs(ius)
ivs=radjacindxs(ivs)
endif
! Initialize ozone jacobian flags to .false. (retain ozone jacobian)
zero_irjaco3_pole = .false.
! These variables are initialized in init_crtm
! isatid = 1 ! index of satellite id
! itime = 2 ! index of analysis relative obs time
! ilon = 3 ! index of grid relative obs location (x)
! ilat = 4 ! index of grid relative obs location (y)
! ilzen_ang = 5 ! index of local (satellite) zenith angle (radians)
! ilazi_ang = 6 ! index of local (satellite) azimuth angle (radians)
! iscan_ang = 7 ! index of scan (look) angle (radians)
! iscan_pos = 8 ! index of integer scan position
! iszen_ang = 9 ! index of solar zenith angle (degrees)
! isazi_ang = 10 ! index of solar azimuth angle (degrees)
! ifrac_sea = 11 ! index of ocean percentage
! ifrac_lnd = 12 ! index of land percentage
! ifrac_ice = 13 ! index of ice percentage
! ifrac_sno = 14 ! index of snow percentage
! its_sea = 15 ! index of ocean temperature
! its_lnd = 16 ! index of land temperature
! its_ice = 17 ! index of ice temperature
! its_sno = 18 ! index of snow temperature
! itsavg = 19 ! index of average temperature
! ivty = 20 ! index of vegetation type
! ivfr = 21 ! index of vegetation fraction
! isty = 22 ! index of soil type
! istp = 23 ! index of soil temperature
! ism = 24 ! index of soil moisture
! isn = 25 ! index of snow depth
! izz = 26 ! index of surface height
! idomsfc = 27 ! index of dominate surface type
! isfcr = 28 ! index of surface roughness
! iff10 = 29 ! index of ten meter wind factor
! ilone = 30 ! index of earth relative longitude (degrees)
! ilate = 31 ! index of earth relative latitude (degrees)
! itref = 34/36 ! index of foundation temperature: Tr
! idtw = 35/37 ! index of diurnal warming: d(Tw) at depth zob
! idtc = 36/38 ! index of sub-layer cooling: d(Tc) at depth zob
! itz_tr = 37/39 ! index of d(Tz)/d(Tr)
! Initialize sensor specific array pointers
! if (goes_img) then
! iclr_sky = 7 ! index of clear sky amount
! elseif (avhrr_navy) then
! isst_navy = 7 ! index of navy sst (K) retrieval
! idata_type = 30 ! index of data type (151=day, 152=night)
! isst_hires = 31 ! index of interpolated hires sst (K)
! elseif (avhrr) then
! iclavr = 32 ! index CLAVR cloud flag with AVHRR data
! isst_hires = 33 ! index of interpolated hires sst (K)
! elseif (seviri) then
! iclr_sky = 7 ! index of clear sky amount
! endif
! Special setup for SST retrieval (output)
if (retrieval.and.init_pass) call setup_sst_retrieval(obstype,dplat(is),mype)
! Special setup for Tz retrieval
if (tzr_qc>0) call setup_tzr_qc(obstype)
! Get version of rad-diag file
call get_radiag ('version',iversion_radiag,istatus)
if(istatus/=0) then
write(6,*)'SETUPRAD: trouble getting version of diag file'
call stop2(999)
endif
! If SSM/I, check for non-use of 85GHz channel, for QC workaround
! set no85GHz true if any 85GHz is not used, and other freq channel is used
no85GHz = .false.
if (ssmi) then
if (iuse_rad(ich(6)) < 1 .or. iuse_rad(ich(7)) < 1 ) then
do j = 1,5
if (iuse_rad(ich(j)) >= 1) then
no85GHz = .true.
cycle
endif
enddo
if (no85GHz .and. mype == 0) write(6,*) &
'SETUPRAD: using no85GHZ workaround for SSM/I ',isis
endif
endif
! Find number of channels written to diag file
if(reduce_diag)then
nchanl_diag=0
do i=1,nchanl
if(iuse_rad(ich(i)) >= 1)then
nchanl_diag=nchanl_diag+1
ich_diag(nchanl_diag)=i
end if
end do
if(mype == mype_diaghdr(is))write(6,*)'SETUPRAD: reduced number of channels ',&
nchanl_diag,' of ',nchanl,' written to diag file '
else
nchanl_diag=nchanl
do i=1,nchanl_diag
ich_diag(i)=i
end do
end if
! Set number of extra pieces of information to write to diagnostic file
! For most satellite sensors there is no extra information. However,
! for GOES Imager data we write additional information.
iextra=0
jextra=0
if (goes_img .or. lwrite_peakwt) then
jextra=nchanl_diag
iextra=1
end if
! If both, iextra=2
if (goes_img .and. lwrite_peakwt) then
iextra=2
end if
lextra = (iextra>0)
! Allocate array to hold channel information for diagnostic file and/or lobsdiagsave option
idiag=ipchan_radiag+npred+3
ioff0 = idiag
if (save_jacobian) then
nnz = nsigradjac
nind = nvarjac
call new(dhx_dx, nnz, nind)
idiag = idiag + size(dhx_dx)
endif
if (lobsdiagsave) idiag=idiag+4*miter+1
allocate(diagbufchan(idiag,nchanl_diag))
allocate(sc_index(nchanl))
sc_index(:) = 0
satinfo_chan: do i=1, nchanl
n = ich(i)
spec_coef: do k=1, sc(1)%n_channels
if ( nuchan(n) == sc(1)%sensor_channel(k)) then
sc_index(i) = k
exit spec_coef
endif
end do spec_coef
end do satinfo_chan
do i=1,nchanl
wavenumber(i)=sc(sensorindex)%wavenumber(sc_index(i))
end do
! If diagnostic file requested, open unit to file and write header.
if (rad_diagsave .and. nchanl_diag > 0) then
if (binary_diag) call init_binary_diag_
if (netcdf_diag) call init_netcdf_diag_
endif
! Load data array for current satellite
read(lunin) data_s,luse,ioid
if (nobskeep>0) then
! write(6,*)'setuprad: nobskeep',nobskeep
call stop2(275)
end if
! PROCESSING OF SATELLITE DATA
! Loop over data in this block
call dtime_setup()
do n = 1,nobs
! Extract analysis relative observation time.
dtime = data_s(itime,n)
call dtime_check(dtime, in_curbin, in_anybin)
if(.not.in_anybin) cycle
if(in_curbin) then
id_qc = igood_qc
if(luse(n))aivals(1,is) = aivals(1,is) + one
! Extract lon and lat.
slons = data_s(ilon,n) ! grid relative longitude
slats = data_s(ilat,n) ! grid relative latitude
cenlon = data_s(ilone,n) ! earth relative longitude (degrees)
cenlat = data_s(ilate,n) ! earth relative latitude (degrees)
! Extract angular information
zasat = data_s(ilzen_ang,n)
cosza = cos(zasat)
zsges=data_s(izz,n)
nadir = nint(data_s(iscan_pos,n))
pangs = data_s(iszen_ang,n)
! Extract warm load temperatures
! wltm1 = data_s(isty,n)
! wltm2 = data_s(istp,n)
! wltm3 = data_s(ism,n)
! If desired recompute 10meter wind factor
if(sfcmod_gfs .or. sfcmod_mm5) then
isli=nint(data_s(idomsfc,n))
call comp_fact10(slats,slons,dtime,data_s(itsavg,n),data_s(isfcr,n), &
isli,mype,data_s(iff10,n))
end if
if(seviri .and. abs(data_s(iszen_ang,n)) > 180.0_r_kind) data_s(iszen_ang,n)=r100
! Set land/sea, snow, ice percentages and flags (no time interpolation)
sea = data_s(ifrac_sea,n) >= 0.99_r_kind
land = data_s(ifrac_lnd,n) >= 0.99_r_kind
ice = data_s(ifrac_ice,n) >= 0.99_r_kind
snow = data_s(ifrac_sno,n) >= 0.99_r_kind
mixed = .not. sea .and. .not. ice .and. &
.not. land .and. .not. snow
eff_area=.false.
if (radmod%lcloud_fwd) then
eff_area=(radmod%cld_sea_only .and. sea) .or. (.not. radmod%cld_sea_only)
end if
iinstr=-1
if(allocated(idnames)) then
if(sea)then
covtype = trim(isis)//':sea'
iinstr=getindex(idnames,trim(covtype))
else if(land)then
covtype = trim(isis)//':land'
iinstr=getindex(idnames,trim(covtype))
else if(ice)then
covtype = trim(isis)//':ice'
iinstr=getindex(idnames,trim(covtype))
else if(snow)then
covtype = trim(isis)//':snow'
iinstr=getindex(idnames,trim(covtype))
else if(mixed)then
covtype = trim(isis)//':mixed'
iinstr=getindex(idnames,trim(covtype))
endif
endif
do jc=1,nchanl
j=ich(jc)
tnoise(jc)=varch(j)
if(sea .and. (varch_sea(j)>zero)) tnoise(jc)=varch_sea(j)
if(land .and. (varch_land(j)>zero)) tnoise(jc)=varch_land(j)
if(ice .and. (varch_ice(j)>zero)) tnoise(jc)=varch_ice(j)
if(snow .and. (varch_snow(j)>zero)) tnoise(jc)=varch_snow(j)
if(mixed .and. (varch_mixed(j)>zero)) tnoise(jc)=varch_mixed(j)
tnoise_save = tnoise(jc)
channel_passive=iuse_rad(j)==-1 .or. iuse_rad(j)==0
if (iuse_rad(j)< -1 .or. (channel_passive .and. &
.not.rad_diagsave)) tnoise(jc)=r1e10
if (passive_bc .and. channel_passive) tnoise(jc)=tnoise_save
if (tnoise(jc) < 1.e4_r_kind) toss = .false.
end do
! Count data of different surface types
if(luse(n))then
if (mixed) then
aivals(5,is) = aivals(5,is) + one
else if (ice .or. snow) then
aivals(4,is) = aivals(4,is) + one
else if (land) then
aivals(3,is) = aivals(3,is) + one
end if
end if
! Set relative weight value
val_obs=one
if(dval_use)then
ixx=nint(data_s(nreal-nstinfo,n))
if (ixx > 0 .and. super_val1(ixx) >= one) then
val_obs=data_s(nreal-nstinfo-1,n)/super_val1(ixx)
endif
endif
! Load channel data into work array.
do i = 1,nchanl
tb_obs(i) = data_s(i+nreal,n)
end do
! Interpolate model fields to observation location, call crtm and create jacobians
! Output both tsim and tsim_clr for allsky
tsim_clr=zero
tcc=zero
if (radmod%lcloud_fwd) then
call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, &
tvp,qvp,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, &
trop5,tzbgr,dtsavg,sfc_speed, &
tsim,emissivity,ptau5,ts,emissivity_k, &
temp,wmix,jacobian,error_status,tsim_clr=tsim_clr,tcc=tcc, &
tcwv=tcwv,hwp_ratio=hwp_ratio,stability=stability)
if(gmi) then
gmi_low_angles(1:3)=data_s(ilzen_ang:iscan_ang,n)
gmi_low_angles(4:5)=data_s(iszen_ang:isazi_ang,n)
data_s(ilzen_ang:iscan_ang, n) = data_s(ilzen_ang2:iscan_ang2, n)
data_s(iszen_ang:isazi_ang, n) = data_s(iszen_ang2:isazi_ang2, n)
call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, &
tvp,qvp,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, &
trop5,tzbgr,dtsavg,sfc_speed, &
tsim2,emissivity2,ptau52,ts2,emissivity_k2, &
temp2,wmix2,jacobian2,error_status,tsim_clr2)
! merge
emissivity(10:13) = emissivity2(10:13)
ts(10:13) = ts2(10:13)
emissivity_k(10:13)= emissivity_k2(10:13)
tsim(10:13) = tsim2(10:13)
wmix(:,10:13) = wmix2(:,10:13)
temp(:,10:13) = temp2(:,10:13)
ptau5(:,10:13) = ptau52(:,10:13)
jacobian(:,10:13) = jacobian2(:,10:13)
! ! output angles for channels 1-9
data_s(ilzen_ang:iscan_ang, n) = gmi_low_angles(1:3)
data_s(iszen_ang:isazi_ang, n) = gmi_low_angles(4:5)
tsim_clr(10:13) = tsim_clr2(10:13)
cosza2 = cos(data_s(ilzen_ang2,n))
endif
else
call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, &
tvp,qvp,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, &
trop5,tzbgr,dtsavg,sfc_speed, &
tsim,emissivity,ptau5,ts,emissivity_k, &
temp,wmix,jacobian,error_status)
if(gmi) then
gmi_low_angles(1:3)=data_s(ilzen_ang:iscan_ang,n)
gmi_low_angles(4:5)=data_s(iszen_ang:isazi_ang,n)
data_s(ilzen_ang:iscan_ang, n) = data_s(ilzen_ang2:iscan_ang2, n)
data_s(iszen_ang:isazi_ang, n) = data_s(iszen_ang2:isazi_ang2, n)
call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, &
tvp,qvp,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, &
trop5,tzbgr,dtsavg,sfc_speed, &
tsim2,emissivity2,ptau52,ts2,emissivity_k2, &
temp2,wmix2,jacobian2,error_status)
! merge
emissivity(10:13) = emissivity2(10:13)
ts(10:13) = ts2(10:13)
emissivity_k(10:13)= emissivity_k2(10:13)
tsim(10:13) = tsim2(10:13)
wmix(:,10:13) = wmix2(:,10:13)
temp(:,10:13) = temp2(:,10:13)
ptau5(:,10:13) = ptau52(:,10:13)
jacobian(:,10:13) = jacobian2(:,10:13)
! ! output angles for channels 1-9
data_s(ilzen_ang:iscan_ang, n) = gmi_low_angles(1:3)
data_s(iszen_ang:isazi_ang, n) = gmi_low_angles(4:5)
cosza2 = cos(data_s(ilzen_ang2,n))
endif
endif
! If the CRTM returns an error flag, do not assimilate any channels for this ob
! and set the QC flag to ifail_crtm_qc.
! We currently go through the rest of the QC steps, ensuring that the diagnostic
! files are populated, but this could be changed if it causes problems.
if (error_status == 0) then
varinv(1:nchanl) = val_obs
else
id_qc(1:nchanl) = ifail_crtm_qc
varinv(1:nchanl) = zero
endif
! For SST retrieval, use interpolated NCEP SST analysis
if (retrieval) then
if( avhrr_navy )then
dtp_avh = data_s(idata_type,n)
sstcu=data_s(isst_hires,n) ! not available, assigned as interpolated sst
sstnv=data_s(isst_navy,n)
elseif ( avhrr) then
if ( pangs <= 89.0_r_kind) then ! day time
dtp_avh = 151.0_r_kind
else
dtp_avh = 152.0_r_kind
endif
sstcu=data_s(isst_hires,n) ! not available, assigned as interpolated sst
sstnv=data_s(isst_hires,n) ! not available, assigned as interpolated sst
endif
tsavg5 = data_s(isst_hires,n)
else
tsavg5=data_s(itsavg,n)
tsavg5=tsavg5+dtsavg
endif
! If using adaptive angle dependent bias correction, update the predicctors
! for this part of bias correction. The AMSUA cloud liquid water algorithm
! uses total angle dependent bias correction for channels 1 and 2
if (adp_anglebc) then
do i=1,nchanl
mm=ich(i)
if (goessndr .or. goes_img .or. ahi .or. seviri .or. ssmi .or. ssmis .or. gmi .or. abi) then
pred(npred,i)=nadir*deg2rad
else
pred(npred,i)=data_s(iscan_ang,n)
end if
do j=2,angord
pred(npred-j+1,i)=pred(npred,i)**j
end do
cbias(nadir,mm)=zero
do j=1,angord