-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathread_obs.F90
1898 lines (1740 loc) · 77.9 KB
/
read_obs.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 read_obsmod
!$$$ module documentation block
! . . . .
! module: read_obsmod extra inquire routine for reading obs
! prgmmr: parrish org: np22 date: 2005-06-06
!
! abstract:
!
! program history log:
! 2009-01-05 todling - add gsi_inquire
! 2015-05-01 Liu Ling - Add call to read_rapidscat
! 2015-08-20 zhu - add flexibility for enabling all-sky and using aerosol info in radiance
! assimilation. Use radiance_obstype_search from radiance_mod.
! 2017-05-12 Y. Wang and X. Wang - add dbz to be read in, POC: xuguang.wang@ou.edu
!
! subroutines included:
! sub gsi_inquire - inquire statement supporting fortran earlier than 2003
! sub read_obs - read, select, reformat obs data
!
! Variable Definitions:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
! set default to private
private
! set subroutines to public
public :: gsi_inquire
public :: read_obs
contains
subroutine gsi_inquire (lbytes,lexist,filename,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: gsi_inquire inquire file presence and size
! prgmmr: todling org: np22 date: 2009-01-05
!
! abstract: Inquire file presence and size; to be used when fortran
! 2003 not available or non-compliant.
!
! program history log:
! 2009-01-05 todling
! 2013-05-14 guo - changed the compiler #ifdef dependency on _size_
! inquire to a more portable way.
!
! input argument list:
! mype - mpi task id
! filename - input filename
!
! output argument list:
! lexist - file presence flag
! lbytes - file size (bytes)
!
! attributes:
! language: f90
! machine: Linux-cluster
!
!$$$ end documentation block
use kinds, only: i_kind,i_llong
implicit none
integer(i_llong),intent( out) :: lbytes
logical ,intent( out) :: lexist
character(len=*),intent(in ) :: filename
integer(i_kind) ,intent(in ) :: mype
character(len=256) command, fname
#if defined(__INTEL_COMPILER) && (__INTEL_COMPILER < 1110)
#define __X_OR_Y_OR_Z_FORTRAN_COMPILERS__
#endif
#ifdef __X_OR_Y_OR_Z_FORTRAN_COMPILERS__
#define _DO_NOT_SUPPORT_SIZE_INQUIRE_
#endif
lbytes=-1 ! in case that _size_ specifier is not supported.
#ifndef _DO_NOT_SUPPORT_SIZE_INQUIRE_
inquire(file=trim(filename),exist=lexist,size=lbytes)
! Note that the value of _size_ is defined by Fortran in "file storage units",
! which is not neccesary in byte units. It is not clear if this code had
! assumed the size value to be in byte units, or in whatever units.
#else
inquire(file=trim(filename),exist=lexist)
#endif
if(lexist) then
! Even with a compiler supporting 'size=' specifier, the size value may
! return -1, if a compiler considers that the size can not be determined.
! In that case, the size may be obtained through a user supported
! mechanism, such as reading the size from a system("wc -c") call result.
if(lbytes<0) then
write(fname,'(2a,i4.4)') 'fsize_',trim(filename),mype
write(command,'(4a)') 'wc -c ', trim(filename),' > ', trim(fname)
call system(command)
open(unit=999,file=trim(fname),form='formatted')
read(999,*) lbytes
close(999)
lexist = lbytes>0_i_llong ! skip this file if lbytes <=0
endif
endif
return
end subroutine gsi_inquire
subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_obs_check inquire file presence and size
! prgmmr: todling org: np22 date: 2010-03-05
!
! abstract: Reset file status depending on whether observation time
! matches analysis time and how offtime_date is set. This
! also checks for consistency in satellite data files and
! known types.
! WARNING: some of it looks inconsistent with long-window 4dvar
!
! program history log:
! 2009-xx-xx derber - originally placed inside inquire
! 2009-01-05 todling - move time/type-check out of inquire
! 2010-09-13 pagowski - add anow bufr and one obs chem
! 2013-01-26 parrish - WCOSS debug compile fails with satid not initialized.
! Set satid=1 at start of subroutine to allow debug compile.
! 2013-02-13 eliu - add ssmis
! 2013-07-01 todling/guo - allow user to bypass this check (old bufr support)
! 2014-10-01 ejones - add gmi and amsr2
! 2015-01-16 ejones - add saphir
! 2016-09-19 guo - properly initialized nread, in case of for quick-return cases.
! 2017-11-16 dutta - adding KOMPSAT5 bufr i.d for reading the data.
! 2019-03-27 h. liu - add abi
!
!
! input argument list:
! lexist - file status
! filename - input filename
! jsatid - satellite id
! dtype - satellite type
!
! output argument list:
! lexist - file status
!
! attributes:
! language: f90
! machine: Linux-cluster
!
!$$$ end documentation block
use kinds, only: i_kind,i_llong,r_kind,r_double
use gsi_4dvar, only: iadatebgn,iadateend
use obsmod, only: offtime_data
use convinfo, only: nconvtype,ictype,ioctype,icuse
use chemmod, only : oneobtest_chem,oneob_type_chem,&
code_pm25_ncbufr,code_pm25_anowbufr,code_pm10_ncbufr,code_pm10_anowbufr
implicit none
logical ,intent(inout) :: lexist
character(len=*),intent(in) :: filename
character(len=*),intent(in) :: jsatid
character(len=*),intent(in) :: dtype
integer(i_kind) ,intent(in) :: minuse
integer(i_kind) ,intent(inout) :: nread
integer(i_kind) :: lnbufr,idate,idate2,iret,kidsat
integer(i_kind) :: ireadsb,ireadmg,kx,nc,said
real(r_double) :: satid,rtype
character(8) subset
logical,parameter:: GMAO_READ=.false.
!======================================
! Loon
integer(i_kind),parameter :: loon_id = 599
!======================================
satid=1 ! debug executable wants default value ???
idate=0
nread=0; if(lexist) nread=1 ! in case of a quick return
#ifdef _SKIP_READ_OBS_CHECK_
return
#endif
if(trim(dtype) == 'tcp' .or. trim(filename) == 'tldplrso')return
if(trim(filename) == 'mitmdat' .or. trim(filename) == 'mxtmdat')return
if(trim(filename) == 'satmar')return
if(trim(dtype) == 'dbz' )return
! Use routine as usual
if(lexist .and. trim(dtype) /= 'tcp' )then
lnbufr = 15
open(lnbufr,file=trim(filename),form='unformatted',status ='unknown')
call openbf(lnbufr,'IN',lnbufr)
call datelen(10)
call readmg(lnbufr,subset,idate,iret)
if(iret == 0)then
! Extract date and check for consistency with analysis date
if (idate<iadatebgn.or.idate>iadateend) then
if(offtime_data) then
write(6,*)'***read_obs_check analysis and data file date differ, but use anyway'
else
write(6,*)'***read_obs_check*** ',&
'incompatable analysis and observation date/time',trim(filename),trim(dtype)
lexist=.false.
end if
write(6,*)'Analysis start :',iadatebgn
write(6,*)'Analysis end :',iadateend
write(6,*)'Observation time:',idate
endif
else
write(6,*)'***read_obs_check*** iret/=0 for reading date for ',trim(filename),dtype,jsatid,iret
lexist=.false.
end if
if(lexist)then
if(jsatid == '')then
kidsat=0
else if(jsatid == 'metop-a')then
kidsat=4
else if(jsatid == 'metop-b')then
kidsat=3
else if(jsatid == 'metop-c')then
kidsat=5
else if(jsatid == 'm08')then
kidsat = 55
else if(jsatid == 'm09')then
kidsat = 56
else if(jsatid == 'm10')then
kidsat = 57
else if(jsatid == 'm11')then
kidsat = 70
else if(jsatid == 'n08')then
kidsat=200
else if(jsatid == 'n09')then
kidsat=201
else if(jsatid == 'n10')then
kidsat=202
else if(jsatid == 'n11')then
kidsat=203
else if(jsatid == 'n12')then
kidsat=204
else if(jsatid == 'n14')then
kidsat=205
else if(jsatid == 'n15')then
kidsat=206
else if(jsatid == 'n16')then
kidsat=207
else if(jsatid == 'n17')then
kidsat=208
else if(jsatid == 'n18')then
kidsat=209
else if(jsatid == 'n19')then
kidsat=223
else if(jsatid == 'npp')then
kidsat=224
else if(jsatid == 'n20')then
kidsat=225
else if(jsatid == 'n21')then
kidsat=226
else if(jsatid == 'f08')then
kidsat=241
else if(jsatid == 'f10')then
kidsat=243
else if(jsatid == 'f11')then
kidsat=244
else if(jsatid == 'f13')then
kidsat=246
else if(jsatid == 'f14')then
kidsat=247
else if(jsatid == 'f15')then
kidsat=248
else if(jsatid == 'f16')then
kidsat=249
else if(jsatid == 'trmm')then
kidsat=282
else if(jsatid == 'f17')then
kidsat=285
else if(jsatid == 'f18')then
kidsat=286
else if(jsatid == 'f19')then
kidsat=287
else if(jsatid == 'g08' .or. jsatid == 'g08_prep')then
kidsat=252
else if(jsatid == 'g09' .or. jsatid == 'g09_prep')then
kidsat=253
else if(jsatid == 'g10' .or. jsatid == 'g10_prep')then
kidsat=254
else if(jsatid == 'g11' .or. jsatid == 'g11_prep')then
kidsat=255
else if(jsatid == 'g12' .or. jsatid == 'g12_prep')then
kidsat=256
else if(jsatid == 'g13' .or. jsatid == 'g13_prep')then
kidsat=257
else if(jsatid == 'g14' .or. jsatid == 'g14_prep')then
kidsat=258
else if(jsatid == 'g15' .or. jsatid == 'g15_prep')then
kidsat=259
else if(jsatid == 'g16' .or. jsatid == 'g16_prep')then
kidsat=270
else if(jsatid == 'g17' .or. jsatid == 'g17_prep')then
kidsat=271
else if(jsatid == 'n05')then
kidsat=705
else if(jsatid == 'n06')then
kidsat=706
else if(jsatid == 'n07')then
kidsat=707
else if(jsatid == 'tirosn')then
kidsat=708
else if ( jsatid == 'terra' ) then
kidsat = 783
else if ( jsatid == 'aqua' ) then
kidsat = 784
else if ( jsatid == 'aura' ) then
kidsat = 785
else if ( jsatid == 'gcom-w1' ) then
kidsat = 122
! Temporary comment gpm out here; discrepancy between SAID in bufr file and
! kidsat.
! else if ( jsatid == 'gpm' ) then
! kidsat = 288
else if ( jsatid == 'meghat' ) then
kidsat = 440
else
kidsat = 0
end if
call closbf(lnbufr)
open(lnbufr,file=trim(filename),form='unformatted',status ='unknown')
call openbf(lnbufr,'IN',lnbufr)
call datelen(10)
if(kidsat /= 0)then
lexist = .false.
satloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(ireadsb(lnbufr)==0)then
call ufbint(lnbufr,satid,1,1,iret,'SAID')
end if
if(nint(satid) == kidsat) then
lexist=.true.
exit satloop
end if
nread = nread + 1
end do satloop
else if(trim(filename) == 'prepbufr')then ! RTod: wired-in filename is not a good idea
lexist = .false.
fileloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
call ufbint(lnbufr,rtype,1,1,iret,'TYP')
kx=nint(rtype)
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit fileloop
end if
end do
end do
nread = nread + 1
end do fileloop
else if(trim(filename) == 'wcpbufr')then
lexist = .false.
file2loop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
call ufbint(lnbufr,rtype,1,1,iret,'TYP')
kx=nint(rtype)
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit file2loop
end if
end do
end do
nread = nread + 1
end do file2loop
else if(trim(filename) == 'gps_ref' .or. trim(filename) == 'gps_bnd')then
lexist = .false.
gpsloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(ireadsb(lnbufr)==0)then
call ufbint(lnbufr,satid,1,1,iret,'SAID')
end if
said=nint(satid)
if(((said > 739) .and.(said < 746)).or.(said == 820) .or. (said == 825) .or. &
(said == 786).or. (said == 4) .or.(said == 3).or. &
( GMAO_READ .and. said == 5) .or. &
(said == 421).or. (said == 440).or.(said == 821)) then
lexist=.true.
exit gpsloop
end if
nread = nread + 1
end do gpsloop
else if(trim(filename) == 'prepbufr_profl')then
lexist = .false.
airploop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
call ufbint(lnbufr,rtype,1,1,iret,'TYP')
kx=nint(rtype)
if (trim(dtype)=='uv') then
if (kx==330 .or. kx==430 .or. kx==530) kx=230
if (kx==331 .or. kx==431 .or. kx==531) kx=231
if (kx==332 .or. kx==432 .or. kx==532) kx=232
if (kx==333 .or. kx==433 .or. kx==533) kx=233
if (kx==334 .or. kx==434 .or. kx==534) kx=234
if (kx==335 .or. kx==435 .or. kx==535) kx=235
if (kx==loon_id) kx=loon_id ! Loon
else
if (kx==330 .or. kx==430 .or. kx==530) kx=130
if (kx==331 .or. kx==431 .or. kx==531) kx=131
if (kx==332 .or. kx==432 .or. kx==532) kx=132
if (kx==333 .or. kx==433 .or. kx==533) kx=133
if (kx==334 .or. kx==434 .or. kx==534) kx=134
if (kx==335 .or. kx==435 .or. kx==535) kx=135
end if
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit airploop
end if
end do
end do
nread = nread + 1
end do airploop
else if(trim(filename) == 'satwndbufr')then
lexist = .false.
loop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
! 5 GOES-R AMVs (NC005030, NC005031, NC005032, NC005034 and NC005039)
! are added as the GOES-R bufr file provide do not contain other winds.
! May not be necessary with the operational satwnd BUFR
if(trim(subset) == 'NC005010' .or. trim(subset) == 'NC005011' .or.&
trim(subset) == 'NC005070' .or. trim(subset) == 'NC005071' .or.&
trim(subset) == 'NC005044' .or. trim(subset) == 'NC005045' .or.&
trim(subset) == 'NC005046' .or. trim(subset) == 'NC005064' .or.&
trim(subset) == 'NC005065' .or. trim(subset) == 'NC005066' .or.&
trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or.&
trim(subset) == 'NC005032' .or. trim(subset) == 'NC005034' .or.&
trim(subset) == 'NC005039') then
lexist = .true.
exit loop
endif
nread = nread + 1
end do loop
else if(trim(filename) == 'oscatbufr')then
lexist = .false.
oscatloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(trim(subset) == 'NC012255') then
lexist = .true.
exit oscatloop
endif
end do oscatloop
else if(trim(filename) == 'rapidscatbufr')then
lexist = .false.
rapidscatloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(trim(subset) == 'NC012255') then
lexist = .true.
exit rapidscatloop
endif
nread = nread + 1
end do rapidscatloop
else if(trim(filename) == 'hdobbufr')then
lexist = .false.
loop_hdob: do while(ireadmg(lnbufr,subset,idate2) >= 0)
if(trim(subset) == 'NC004015') then
lexist = .true.
exit loop_hdob
endif
nread = nread + 1
end do loop_hdob
else if(trim(dtype) == 'pm2_5')then
if (oneobtest_chem .and. oneob_type_chem=='pm2_5') then
lexist=.true.
else
lexist = .false.
fileloopanow_pm2_5:do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
if (subset == 'ANOWPM') then
call ufbint(lnbufr,rtype,1,1,iret,'TYP')
kx=nint(rtype)
else if ( (subset == 'NC008031') .or. &
(subset == 'NC008032' ) ) then
call ufbint(lnbufr,rtype,1,1,iret,'TYPO')
kx=nint(rtype)
if (kx/=code_pm25_ncbufr) then
cycle
else
kx=code_pm25_anowbufr
endif
else
cycle
endif
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. &
kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit fileloopanow_pm2_5
end if
end do
end do
nread = nread + 1
enddo fileloopanow_pm2_5
endif
if (lexist) then
write(6,*)'found pm2_5 in anow bufr'
else
write(6,*)'did not find pm2_5 in anow bufr'
endif
else if(trim(dtype) == 'pm10')then
lexist = .false.
fileloopanow_pm10:do while(ireadmg(lnbufr,subset,idate2) >= 0)
do while(ireadsb(lnbufr)>=0)
if (subset == 'NC008033') then
call ufbint(lnbufr,rtype,1,1,iret,'TYPO')
kx=nint(rtype)
IF (kx/=code_pm10_ncbufr) then
cycle
else
kx=code_pm10_anowbufr
endif
else
cycle
endif
do nc=1,nconvtype
if(trim(ioctype(nc)) == trim(dtype) .and. &
kx == ictype(nc) .and. icuse(nc) > minuse)then
lexist = .true.
exit fileloopanow_pm10
end if
end do
end do
nread = nread + 1
enddo fileloopanow_pm10
if (lexist) then
write(6,*)'found pm10 in anow bufr'
else
write(6,*)'did not find pm10 in anow bufr'
endif
end if
end if
call closbf(lnbufr)
end if
if(lexist)then
write(6,*)'read_obs_check: bufr file date is ',idate,trim(filename),' ',dtype,jsatid
else
write(6,*)'read_obs_check: bufr file ',dtype,jsatid,' not available ',trim(filename)
end if
return
end subroutine read_obs_check
subroutine read_obs(ndata,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_obs read, select, reformat obs data
! prgmmr: parrish org: np22 date: 1990-10-07
!
! abstract: This routine is a driver for routines which read different
! types of observational data.
!
! program history log:
! 1990-10-07 parrish
! 1998-05-15 weiyu yang
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2004-06-16 treadon - update documentation
! 2004-07-23 derber - modify to include conventional sst
! 2004-07-29 treadon - add only to module use, add intent in/out
! 2005-01-20 okamoto - add calling read_ssmi
! 2005-06-14 wu - add OMI oz
! 2005-07-06 derber - add mhs, hirs/4 and ears data
! 2005-08-02 derber - modify to use convinfo file
! 2005-09-08 derber - modify to use input group time window
! 2005-09-20 xu & pawlak - modify calling read_ssmis and read_amsre
! 2005-09-28 derber - modify to simplify obs handling
! 2005-10-17 derber - pass obs_load1 into read_amsre and read_ssmis
! 2005-10-18 treadon - remove obs_load and obs_load1
! 2005-10-20 kazumori - modify to read real AMSR-E data
! 2005-11-28 derber move determination of which ob data sets to read inside read_obs
! 2005-11-14 li, xu - modify sst obs read and add avhrr gac 1b obs read
! 2006-01-25 treadon - remove read_ieeetovs
! 2006-02-01 parrish - add getsfc and destroy_sfc for full surface fields
! 2006-02-03 derber - modify for new obs control and obs count
! 2005-02-03 treadon - gather guess 3d pressure to full grid array
! 2006-02-08 derber - modify to use new convinfo module
! 2006-02-01 liu - add ssu
! 2006-02-24 derber - modify to take advantage of convinfo module
! 2006-03-07 derber - consolidate processing of 1x1 and 5x5 goessndr
! 2006-04-20 kistler - moved conv_read to gsisub
! 2006-05-25 treadon - rename goesimg and goes_img and pcp_ssm/i as
! pcp_ssmi to make consistent with other obstype
! 2006-09-20 treadon - add mpi_io for select data file (obstype)s
! 2006-10-12 treadon - remove tendsflag check for pcp data (now in gsimain)
! 2007-02-21 sienkiewicz - bring in changes for MLS ozone
! 2007-03-15 su - add reading conventional error table option
! 2007-06-05 treadon - restructure mpi_querybf section to improve efficiency
! 2007-10-03 todling - skip most of this in 4dvar inner loop
! 2008-03-28 wu - move random seed for perturb_obs from setuprhsall
! 2008-04-18 safford - rm unused vars and uses
! 2008-05-01 h.liu - add gome ozone
! 2008-06-20 derber - move destroy_sfc to this routine
! 2008-09-08 lueken - merged ed''s cahnges into q1fy09 code
! 2008-12-30 todling - handle inquire for diff versions of fortran
! 2009-01-05 todling - need tendency alloc in observer mode
! 2009-01-23 todling - echo surface state info
! 2009-03-18 meunier - add a if statement to read lagrangian data
! 2009-08-19 guo - moved destroy_sfc_grid() to observer_finalize().
! 2009-12-20 gayno - modify argument lists so that fov-based surface
! calculation may be used.
! 2010-03-29 hu - add code to read in cloud observations including:
! prepbufr (metar, nesdis cloud product)
! radar reflectivity, lightning, NASA LaRC cloud
! 2010-04-01 treadon - move strip and reorder to gridmod
! 2010-04-08 hliu - add seviri
!
! 2010-04-05 huang - add aero and modis for reading modis aod from satellite
! currently read BUFR file only
! 2010-04-22 tangborn - read for carbon monoxide
! 2010-08-23 tong - add calcuation of hgtl_full used in read_radar.f90
! 2011-04-02 li - add nst_gsi, getnst and destroy_nst
! 2011-05-20 mccarty - add cris/atms handling
! 2011-05-26 todling - add call to create_nst
! 2012-03-07 akella - set NST variables; init & destroy in guess_grids
! [see create_sfc_grids & destroy_sfc_grids]
! 2012-05-16 todling - protect call to prebufr when NST is on (should be mods)
! 2012-09-08 wargan - add OMI with efficiency factors
! 2013-01-26 parrish - WCOSS debug compile fails--extra arguments in call read_aerosol.
! Commented out extra line of arguments not used.
! 2013-02-13 eliu - turn off parallel I/O for SSMIS (due to the need to
! do spatial averaging for noise reduction)
! 2013-06-01 zhu - add mype_airobst to handle aircraft temperature bias correction
! 2013-08-08 s.liu - add read NASA_LaRC_cloud product
! 2013-10-25 todling - reposition ltosi and others to commvars
! 2014-01-01 xli - add option to read NSST marine BUFR data file nsstbufr (on the
! top of prepbufr and modsbufr)
! 2014-02-03 guo - Hid processing (read) of non-EMC ozone obstypes
! through module m_extOzone, separated from read_ozone.
! - Added some -do- and -if- construct names, for easier
! understanding of the code.
! 2014-06-19 carley/zhu - Add tcamt and lcbas
! 2014-11-12 carley - Add call to read goes imager sky cover data for tcamt
! 2014-12-03 derber - modify for no radiance cases and read processor for
! surface fields
! 2015-01-16 ejones - added saphir, gmi, and amsr2 handling
! 2015-03-23 zaizhong ma - add Himawari-8 ahi
! 2015-05-30 li - modify for no radiance cases but sst (nsstbufr) and read processor for
! surface fields (use_sfc = .true. for data type of sst),
! to use deter_sfc in read_nsstbufr.f90)
! 2015-07-10 pondeca - add cloud ceiling height (cldch)
! 2015-08-12 pondeca - add capability to read min/maxT obs from ascii file
! 2015-08-20 zhu - use centralized radiance info from radiance_mod to specify
! all-sky and aerosol usages in radiance assimilation
! 2015-09-04 J. Jung - Added mods for CrIS full spectral resolution (FSR)
! 2015-10-19 lippi - Added logic to read l2rw or rw radial winds.
! 2016-03-02 s.liu/carley - remove use_reflectivity and use i_gsdcldanal_type
! 2016-04-28 J. Jung - added logic for RARS and direct broadcast data from NESDIS/UW.
! 2016-05-05 pondeca - add 10-m u-wind and v-wind (uwnd10m, vwnd10m)
! 2016-09-19 Guo - replaced open(obs_input_common) with "call unformatted_open(obs_input_common)"
! 2017-05-12 Y. Wang and X. Wang - add multi-interface to read in dBZ (nc) and radial velocity (ascii)
! 2016-12-14 lippi - Fixed bug of using observations twice when both
! l2rwbufr and radarbufr are in the OBS_INPUT table.
! Changed the dsis entries for l2rwbufr and radarbufr to
! l2rw and l3rw respectively. Also make use of nml
! option vadwnd_l2rw_qc.
! 2017-08-31 Li - move gsi_nstcoupler_init & gsi_nstcoupler_set to getsfc in sathin.F90
! - move gsi_nstcoupler_final from create_sfc_grids to here
! 2017-12-05 Wargan - added OMPS ozone
! 2018-01-23 Apodaca - add GOES/GLM lightning data
! 2018-07-09 Todlng - move gsi_nstcoupler_final to destroy_sfc (consistency)
! 2019-01-15 Li - add to handle mbuoyb
! 2019-03-27 h. liu - add abi
!
!
! input argument list:
! mype - mpi task id
!
! output argument list:
! ndata(*,1)- number of prefiles retained for further processing
! ndata(*,2)- number of observations read
! ndata(*,3)- number of observations keep after read
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind,i_llong,r_double
use gridmod, only: nlon,nlat,nsig,iglobal,ijn,itotsub,lat1,lon1,&
displs_g,strip,reorder
use general_commvars_mod, only: ltosi,ltosj
use obsmod, only: iadate,ndat,time_window,dplat,dsfcalc,dfile,dthin, &
dtype,dval,dmesh,obsfile_all,ref_obs,nprof_gps,dsis,ditype,&
perturb_obs,lobserver,lread_obs_save,obs_input_common, &
reduce_diag,nobs_sub,dval_use
use gsi_nstcouplermod, only: nst_gsi
! use gsi_nstcouplermod, only: gsi_nstcoupler_set
use qcmod, only: njqc,vadwnd_l2rw_qc
use gsi_4dvar, only: l4dvar
use satthin, only: super_val,super_val1,superp,makegvals,getsfc,destroy_sfc
use mpimod, only: ierror,mpi_comm_world,mpi_sum,mpi_rtype,mpi_integer,npe,&
setcomm
use constants, only: one,zero
use converr, only: converr_read
use converr_ps, only: converr_ps_read
use converr_q, only: converr_q_read
use converr_t, only: converr_t_read
use converr_uv, only: converr_uv_read
use converr_pw, only: converr_pw_read
use convb_ps,only: convb_ps_read
use convb_q,only:convb_q_read
use convb_t,only:convb_t_read
use convb_uv,only:convb_uv_read
use guess_grids, only: ges_prsl,geop_hgtl,ntguessig
use radinfo, only: nusis,iuse_rad,jpch_rad,diag_rad
use insitu_info, only: mbuoy_info,mbuoyb_info,read_ship_info
use aeroinfo, only: nusis_aero,iuse_aero,jpch_aero,diag_aero
use ozinfo, only: nusis_oz,iuse_oz,jpch_oz,diag_ozone
use pcpinfo, only: npcptype,nupcp,iusep,diag_pcp
use convinfo, only: nconvtype,ioctype,icuse,diag_conv,ithin_conv
use chemmod, only : oneobtest_chem,oneob_type_chem,oneobschem
use lightinfo, only: nlighttype,iuse_light,diag_light
use aircraftinfo, only: aircraft_t_bc,aircraft_t_bc_pof,aircraft_t_bc_ext,mype_airobst
use gsi_io, only: mype_io
use rapidrefresh_cldsurf_mod, only: i_gsdcldanal_type
use radiance_mod, only: rad_obs_type,radiance_obstype_search
use m_extOzone, only: is_extOzone
use m_extOzone, only: extOzone_read
use mpeu_util, only: warn
use gsi_unformatted, only: unformatted_open
use mrmsmod,only: l_mrms_sparse_netcdf
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: mype
integer(i_kind),dimension(ndat,3),intent( out) :: ndata
! Declare local parameters
integer(i_llong),parameter:: lenbuf=8388608_i_llong ! lenbuf=8*1024*1024
! Declare local variables
logical :: lexist,ssmis,amsre,sndr,hirs,avhrr,lexistears,lexistdb,use_prsl_full,use_hgtl_full
logical :: use_sfc,nuse,use_prsl_full_proc,use_hgtl_full_proc,seviri,mls,abi
logical,dimension(ndat):: belong,parallel_read,ears_possible,db_possible
logical :: modis,use_sfc_any
logical :: acft_profl_file
character(10):: obstype,platid
character(22):: string
character(120):: infile
character(20):: sis
integer(i_kind) i,j,k,ii,nmind,lunout,isfcalc,ithinx,ithin,nread,npuse,nouse
integer(i_kind) nprof_gps1,npem1,krsize,len4file,npemax,ilarge,nlarge,npestart
integer(i_llong) :: lenbytes
integer(i_kind):: npetot,npeextra,mmdat,nodata
integer(i_kind):: iworld,iworld_group,next_mype,mm1,iix
integer(i_kind):: mype_root
integer(i_kind):: minuse,lunsave,maxproc,minproc
integer(i_kind),dimension(ndat):: npe_sub,npe_sub3,mpi_comm_sub,mype_root_sub,npe_order
integer(i_kind),dimension(ndat):: ntasks1,ntasks
integer(i_kind),dimension(ndat):: read_rec1,read_rec
integer(i_kind),dimension(ndat):: read_ears_rec1,read_ears_rec
integer(i_kind),dimension(ndat):: read_db_rec1,read_db_rec
integer(i_kind),dimension(ndat,3):: ndata1
integer(i_kind),dimension(npe,ndat):: mype_work,nobs_sub1
integer(i_kind),dimension(npe,ndat):: mype_sub
integer(i_kind),allocatable,dimension(:):: nrnd
integer(i_kind):: nmls_type,mype_io_sfc
integer(i_kind):: iread,ipuse,iouse
real(r_kind) gstime,val_dat,rmesh,twind,rseed
real(r_kind),allocatable,dimension(:) :: prslsm,hgtlsm,work1
real(r_kind),allocatable,dimension(:,:,:):: prsl_full,hgtl_full
type(rad_obs_type) :: radmod
data lunout / 81 /
data lunsave / 82 /
!*****************************************************************************
! Set analysis time and allocate/initialize arrays and variables
call w3fs21(iadate,nmind)
gstime=real(nmind,r_kind)
allocate(nobs_sub(npe,ndat))
nobs_sub = 0
nobs_sub1 = 0
call makegvals
do ii=1,ndat
ndata1(ii,1)=0
ndata1(ii,2)=0
ndata1(ii,3)=0
ntasks1(ii) =0
parallel_read=.false.
end do
npem1=npe-1
nprof_gps1=0
!obs error tables are read here!
if(njqc) then
call converr_ps_read(mype)
call converr_q_read(mype)
call converr_t_read(mype)
call converr_uv_read(mype)
call converr_pw_read(mype)
call convb_ps_read(mype)
call convb_q_read(mype)
call convb_t_read(mype)
call convb_uv_read(mype)
else
call converr_read(mype)
endif
! Optionally set random seed to perturb observations
if (perturb_obs) then
rseed=iadate(4)+iadate(3)*100+iadate(2)*10000+iadate(1)*1000000+mype
call random_seed(size=krsize)
allocate(nrnd(krsize))
do i=1,krsize
nrnd(i)=rseed
end do
call random_seed(put=nrnd)
deallocate(nrnd)
endif
! Set data class and number of reader tasks. Set logical flag to indicate
! type type of GPS data (if present)
ii=0
ref_obs = .false. !.false. = assimilate GPS bending angle
ears_possible = .false.
db_possible = .false.
nmls_type=0
read_rec1 = 0
read_ears_rec1=0
read_db_rec1=0
do i=1,ndat
obstype=dtype(i) ! obstype - observation types to process
amsre= index(obstype,'amsre') /= 0
ssmis= index(obstype,'ssmis') /= 0
sndr = index(obstype,'sndr') /= 0
hirs = index(obstype,'hirs') /= 0
avhrr = index(obstype,'avhrr') /= 0
modis = index(obstype,'modis') /= 0
seviri = index(obstype,'seviri') /= 0
abi = index(obstype,'abi') /= 0
mls = index(obstype,'mls') /= 0
if(obstype == 'mls20' ) nmls_type=nmls_type+1
if(obstype == 'mls22' ) nmls_type=nmls_type+1
if(obstype == 'mls30' ) nmls_type=nmls_type+1
if(nmls_type>1) then
write(6,*) '******ERROR***********: there is more than one MLS data type, not allowed, please check'
call stop2(339)
end if
if (obstype == 't' .or. obstype == 'uv' .or. &
obstype == 'q' .or. obstype == 'ps' .or. &
obstype == 'pw' .or. obstype == 'spd'.or. &
obstype == 'sst'.or. &
obstype == 'tcp'.or. obstype == "lag".or. &
obstype == 'dw' .or. obstype == 'rw' .or. &
obstype == 'mta_cld' .or. obstype == 'gos_ctp' .or. &
obstype == 'rad_ref' .or. obstype=='lghtn' .or. &
obstype == 'larccld' .or. obstype == 'pm2_5' .or. obstype == 'pm10' .or. &
obstype == 'gust' .or. obstype=='vis' .or. &
obstype == 'pblh' .or. obstype=='wspd10m' .or. &
obstype == 'td2m' .or. obstype=='mxtm' .or. &
obstype == 'mitm' .or. obstype=='pmsl' .or. &
obstype == 'howv' .or. obstype=='tcamt' .or. &
obstype=='lcbas' .or. obstype=='cldch' .or. obstype == 'larcglb' .or. &
obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' ) then
ditype(i) = 'conv'
else if (obstype == 'swcp' .or. obstype == 'lwcp') then
ditype(i) = 'wcp'
else if( hirs .or. sndr .or. seviri .or. abi .or. &
obstype == 'airs' .or. obstype == 'amsua' .or. &
obstype == 'msu' .or. obstype == 'iasi' .or. &
obstype == 'amsub' .or. obstype == 'mhs' .or. &
obstype == 'hsb' .or. obstype == 'goes_img' .or. &
obstype == 'ahi' .or. avhrr .or. &
amsre .or. ssmis .or. obstype == 'ssmi' .or. &
obstype == 'ssu' .or. obstype == 'atms' .or. &
obstype == 'cris' .or. obstype == 'cris-fsr' .or. &
obstype == 'amsr2' .or. &
obstype == 'gmi' .or. obstype == 'saphir' ) then
ditype(i) = 'rad'
else if (is_extOzone(dfile(i),obstype,dplat(i))) then
ditype(i) = 'ozone'
else if (obstype == 'sbuv2' &
.or. obstype == 'omi' &
.or. obstype == 'ompstc8' &
.or. obstype == 'ompsnp' &
.or. obstype == 'gome' &
.or. index(obstype, 'omps') /= 0 &
.or. mls &
) then
ditype(i) = 'ozone'
else if (obstype == 'mopitt') then
ditype(i) = 'co'
else if (index(obstype,'pcp')/=0 )then
ditype(i) = 'pcp'
else if (obstype == 'gps_ref' .or. obstype == 'gps_bnd') then
ditype(i) = 'gps'
else if ( index(obstype,'aod') /= 0 ) then
ditype(i) = 'aero'
else if (obstype == 'goes_glm') then
ditype(i) = 'light'
else
write(6,*)'READ_OBS: ***ERROR*** - unknown ob type ',trim(obstype)
end if
! Set data class and number of reader tasks. Set logical flag to indicate
! type type of GPS data (if present)
if (index(dtype(i),'gps_ref') /= 0) ref_obs = .true.
! Check info files to see if data is read.
nuse=.false.
minuse=-1
if(ditype(i) == 'conv')then
if(diag_conv .and. .not. reduce_diag)minuse=-2
do j=1,nconvtype
if(trim(dtype(i)) == trim(ioctype(j)) .and. icuse(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'rad')then
if(diag_rad .and. .not. reduce_diag)minuse=-2
do j=1,jpch_rad
if(trim(dsis(i)) == trim(nusis(j)) .and. iuse_rad(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'ozone')then
if(diag_ozone .and. .not. reduce_diag)minuse=-2
do j=1,jpch_oz
if(trim(dsis(i)) == trim(nusis_oz(j)) .and. iuse_oz(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'pcp')then
if(diag_pcp .and. .not. reduce_diag)minuse=-2
do j=1,npcptype
if(trim(dsis(i)) == trim(nupcp(j)) .and. iusep(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'aero')then
if(diag_aero .and. .not. reduce_diag)minuse=-2
do j=1,jpch_aero
if(trim(dsis(i)) == trim(nusis_aero(j)) .and. iuse_aero(j) > minuse)nuse=.true.
end do
else if(ditype(i) == 'light')then
if(diag_light .and. .not. reduce_diag)minuse=-2
do j=1,nlighttype
if(iuse_light(j) > minuse) nuse=.true.
end do
else
nuse=.true.
end if
if(nuse)then
! Control parallel read for each ob type (currently just rad obs).
! To remove parallel read comment out line.
ithin=dthin(i)
if(ithin > 0 )then
if(dmesh(ithin) > one)then
if(hirs)then
parallel_read(i)= .true.
else if(obstype == 'amsua')then
parallel_read(i)= .true.
else if(obstype == 'airs' )then
parallel_read(i)= .true.
else if(obstype == 'iasi')then
parallel_read(i)= .true.
else if(obstype == 'amsub')then
parallel_read(i)= .true.
else if(obstype == 'mhs' )then
parallel_read(i)= .true.
else if(sndr )then
parallel_read(i)= .true.
! N.B. ATMS must be run on one processor for the filtering code to work.
else if(obstype == 'atms')then
! parallel_read(i)= .true.
else if(ssmis)then
! parallel_read(i)= .true.
else if(seviri)then
parallel_read(i)= .true.
else if(abi)then
parallel_read(i)= .true.
else if(obstype == 'cris' .or. obstype == 'cris-fsr')then
parallel_read(i)= .true.
else if(avhrr)then
parallel_read(i)= .true.
else if(amsre)then
parallel_read(i)= .true.
else if(obstype == 'goes_img' )then
parallel_read(i)= .true.
else if(obstype == 'ahi' )then