forked from DOI-USGS/volcano-ash3d-metreader
-
Notifications
You must be signed in to change notification settings - Fork 4
/
MetReader_ASCII.f90
3846 lines (3658 loc) · 268 KB
/
MetReader_ASCII.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
! Here is a script for downloading the Radiosonde data in Alaska
!#!/bin/bash
!
!StationNum=("70414" "70316" "70326" "70350" "70273")
!StationCode=("PASY" "PACD" "PAKN" "PADQ" "PANC")
!StationName=("Shemya Afb" "Cold Bay" "King Salmon" "Kodiak" "Anchorage")
!#yearmonthday=`date -u +%Y%m%d` #current year, month &
!day (e.g. 20110119)
!Y=`date -u +%Y`
!M=`date -u +%m`
!DD=`date -u +%d`
!Y=2015
!M=02
!DD=24
!HH=(00 12)
!
!cd /data/windfiles/MetProfiles
!
!for (( si=0;si<=4;si++))
!do
! #TYPE=TEXT%3AUNMERGED
! echo
!"http://weather.uwyo.edu/cgi-bin/sounding?region=naconf&TYPE=TEXT%3ALIST&YEAR=${Y}&MONTH=${M}&FROM=${DD}${HH[$1]}&TO=${DD}${HH[$1]}&STNM=${StationNum[si]}"
!> tmp.lnk
! wget -i tmp.lnk -O out.html
! cp out.html RAW/${StationCode[si]}_${Y}${M}${DD}${HH[$1]}_raw.dat
!
! sed '1,10d' out.html > headless.dat
! sed '/PRE/,$d' headless.dat > ${StationCode[si]}_${Y}${M}${DD}${HH[$1]}.dat
! rm tmp.lnk headless.dat out.html
!
!done
!
!PRES: Atmospheric Pressure [hPa]
!HGHT: Geopotential Height [meter]
!TEMP: Temperature [Celsius]
!DWPT: Dewpoint Temperature [Celsius]
!RELH: Relative Humidity [%]
!MIXR: Mixing Ratio [gram/kilogram]
!DRCT: Wind Direction [degrees true]
!SKNT: Wind Speed [knot]
!THTA: Potential Temperature [kelvin]
!THTE: Equivalent Potential Temperature [kelvin]
!THTV: Virtual Potential Temperature [kelvin]
!##############################################################################
!
! MR_Read_Met_DimVars_ASCII_1d
!
! Called once from MR_Read_Met_DimVars
!
! This subroutine needs to determine the full Met grid (i.e. determine the
! pressure grid and the spatial grid). It also needs to determine which
! variables are available. Generally, the actual region used (e.g. the submet)
! is determined through MR_Initialize_Met_Grids after the comp grid is
! defined. Values are normally only read on this submet grid. For 1d
! ASCII cases, all sonde data are loaded into memory in this subroutine.
! MR_Initialize_Met_Grids will be used to define the mapping of sonde data
! to the comp grid.
!
! This subroutine also allocates and populates the time arrays for the sonde
! data, something usually done by MR_Read_Met_Times_[netcdf,grib,etc]
! Note: Every sonde profile will have its own pressure levels, however,
! the 1d array, p_fullmet_sp, is allocated as the intersection of all
! the individual pressure profiles so that programs can have data on the
! metP grid
!
! Allocated the dummy arrays for storing met data on met and computational
! grids: MR_dum2d_met(nx_submet,ny_submet)
! MR_dum3d_metP(nx_submet,ny_submet,np_fullmet)
! MR_dum3d_metH(nx_submet,ny_submet,nz_comp)
! MR_dum2d_comp(nx_comp,ny_comp)
! MR_dum3d_compH(nx_comp,ny_comp,nz_comp)
! MR_windfile_starthour(MR_Snd_nt_fullmet)
! MR_windfile_stephour(MR_Snd_nt_fullmet,1)
!
!##############################################################################
subroutine MR_Read_Met_DimVars_ASCII_1d
use MetReader, only : &
MR_nio,VB,outlog,errlog,verbosity_error,verbosity_info,verbosity_production,verbosity_debug1,&
x_fullmet_sp,y_fullmet_sp,p_fullmet_sp,Snd_idx,MR_nSnd_Locs,MR_iwindfiles,MR_windfile_stephour,&
MR_windfiles_nt_fullmet,MR_windfile_starthour,Met_dim_IsAvailable,Met_var_IsAvailable,&
MR_SndVars_metP,MR_SndVarsID,MR_Snd_np_fullmet,np_fullmet,z_approx,MR_windfiles,&
Snd_Have_PT,nt_fullmet,MR_BaseYear,MR_useLeap,MR_Snd_nt_fullmet,MR_Snd_nvars,Snd_Have_Coord,&
MR_Max_geoH_metP_predicted,MR_iwind,MR_iwindformat,IsGlobal_MetGrid,IsLatLon_MetGrid,IsRegular_MetGrid,&
Met_iprojflag,Met_k0,Met_lam0,Met_lam1,Met_lam2,Met_phi0,Met_phi1,Met_phi2,Met_Re,MR_EPS_SMALL,&
MR_Comp_StartYear,MR_Comp_StartMonth,MR_nstat,MR_dx_met,MR_dy_met,&
MR_Z_US_StdAtm,&
MR_Temp_US_StdAtm,&
MR_Pres_US_StdAtm, &
MR_FileIO_Error_Handler
use projection, only : &
PJ_iprojflag,PJ_k0,PJ_lam0,PJ_lam1,PJ_lam2,PJ_phi0,PJ_phi1,PJ_phi2,PJ_Re,&
PJ_Set_Proj_Params
implicit none
integer, parameter :: sp = 4 ! single precision
integer, parameter :: dp = 8 ! double precision
real(kind=dp), parameter :: PI = 3.141592653589793_dp
real(kind=dp), parameter :: DEG2RAD = 1.7453292519943295e-2_dp
real(kind=sp), parameter :: KNOTS2MS = 0.514444444_sp
integer, parameter :: MAX_ROWS = 300 ! maximum number of row of data
integer, parameter :: fid = 120
integer :: iostatus
integer :: ioerr
character(len=120) :: iomessage = ""
integer :: istr1,istr2,istr3
integer :: ic,il,iil,iv
integer :: nlev,ulev
integer :: iw_idx
integer :: iloc, itime
integer :: p_lidx, p_tidx
real(kind=sp) :: p_maxtop, p_top
integer :: k
real(kind=sp),dimension(:),allocatable :: WindVelocity
real(kind=sp),dimension(:),allocatable :: WindDirection
real(kind=sp),dimension(16) :: pres_Snd_tmp
real(kind=sp) :: WindTime
integer :: iw,iws
real(kind=sp) :: rvalue1,rvalue2,rvalue3
real(kind=sp),dimension(10) :: rvalues
integer :: ivalue1,ivalue2,ivalue3,ivalue4,ivalue5
real(kind=sp) :: rvalue1_o,rvalue3_o
integer :: ivalue2_o,ivalue4_o,ivalue5_o
integer :: ilatlonflag
character(len=80) :: linebuffer080
character(len=80) :: linebuffer080_2
character(len=80) :: linebuffer080_3
character(len=80) :: linebuffer080_4
character(len=7) :: field_str
character(len=6),dimension(53) :: GTSstr
character(len=6) :: dumstr1,dumstr2,dumstr3,dumstr4
real(kind=sp) :: WindSuppl
integer :: dum_int
integer,dimension(:),allocatable :: H_tmp
integer,dimension(:),allocatable :: T_tmp
integer :: il_bot,il_top,ill
real(kind=sp) :: vbot,vtop,zbot,ztop,zhere
integer :: scl_idx
real(kind=sp),dimension(7) :: scl_m,scl_a
real(kind=sp) :: SurfPres
real(kind=sp) :: SurfTemp
real(kind=sp) :: SurfDewPoint
real(kind=sp) :: SurfWindDir
real(kind=sp) :: SurfWindSpeed
integer :: SurfTemp_int
logical :: In_hPa = .true.
integer :: indx1,indx2,indx3
integer :: ncols
logical :: IsWindDirectSpeed
logical :: IsCustVarOrder
integer,dimension(:),allocatable :: SndColReadOrder
character(len=3),dimension(0:50) :: MR_SndVarsName
integer :: Stat_ID
integer :: Stat_idx
real(kind=sp) :: Stat_elev
integer :: idx,idx2
real(kind=sp) :: fac
integer :: substr_pos1
character(len=8) :: date
character(len=10) :: time2
character(len=5) :: zone
integer :: values(8)
logical :: IsGTS = .false.
logical :: IsRUCNOAA = .false.
logical :: HasTTCC = .false.
integer :: DayOfMonth, SndHour
integer :: io ! Index for output streams
INTERFACE
subroutine MR_Load_Radiosonde_Station_Data
end subroutine MR_Load_Radiosonde_Station_Data
subroutine MR_Get_Radiosonde_Station_Coord(StatID,StatIdx,Stat_lon,Stat_lat,Stat_elv)
integer, parameter :: sp = 4 ! single precision
integer, intent(in) :: StatID
integer, intent(out) :: StatIdx
real(kind=sp),intent(out) :: Stat_lon,Stat_lat,Stat_elv
end subroutine MR_Get_Radiosonde_Station_Coord
real(kind=8) function HS_hours_since_baseyear(iyear,imonth,iday,hours,byear,useLeaps)
integer ,intent(in) :: iyear
integer ,intent(in) :: imonth
integer ,intent(in) :: iday
real(kind=8),intent(in) :: hours
integer ,intent(in) :: byear
logical ,intent(in) :: useLeaps
end function HS_hours_since_baseyear
END INTERFACE
do io=1,MR_nio;if(VB(io).le.verbosity_production)then
write(outlog(io),*)"-----------------------------------------------------------------------"
write(outlog(io),*)"---------- MR_Read_Met_DimVars_ASCII_1d ----------"
write(outlog(io),*)"-----------------------------------------------------------------------"
endif;enddo
!------------------------------------------------------------------------------
! MR_iwind.eq.1
!--------------------------------------------------------------------
! MR_iwindformat.eq.1 (custom 1-d profile)
! L1 string header line
! L2 time(hr) nlev [ncol] [ivar(ncol)]
! L3 x/lon y/lat [LLflag] []
! data
! if no optional params, then
! 3 col -> alt(m) windspeed(m/s) winddir(E of N wind FROM)
! 5 col -> alt(m) windspeed(m/s) winddir(E of N wind FROM) pres(hPa) T(degC)
! ncol -> variables can be in any order but must include:
! 1) either pres or alt or both (if one is given, the other is filled from stdAtm)
! 2) either vx/vy (vector componets) or speed/dir (with 'dir' as wind from)
! If T is not provided, it is auto-filled from StdAtm
!-------------
!Input wind file for Ashfall tests
!2 41 #wind time, number of levels
!-122.18 46.20 #Longitude, latitude of wind sounding (fictitious)
!0 10.000 90.00 101300. 15.0 z, windspeed, direction, p, T
!1000 10.000 90.00 89846. 8.5
!2000 10.000 90.00 79464. 2.0
!--------------------------------------------------------------------
! MR_iwindformat.eq.2 (radiosonde format as downloaded from "http://weather.uwyo.edu)
! This can either be in the Raw or WMO/GTS format, or the 'list' format.
! Example WMO/GTS
! TTAA 56001 72694 99009 22868 29006 00137 20667 28507 92801 13662 25507 85505
!06857 21003 70077 01471 23023 50573 12971 24538 40739 25780 25040 30939 43564
!23045 25060 52563 23558 20200 61761 22548 15380 565// 24035 10637 567// 22531
!88185 62963 23051 77221 23563 40311 31313 58208 82326 51515 10164 00012 10194
!26506 22517
! TTBB 56008 72694 00009 22868 11999 20467 22979 18466 33806 02640 44802 02444
!55797 03459 66793 03658 77784 03256 88776 03456 99758 02057 11751 03070 22739
!02890 33717 01481 44713 01672 55702 01273 66698 01468 77665 00760 88660 00563
!99651 01161 11646 01165 22636 01170 33629 01365 44610 02363 55599 02965 66550
!07169 77523 10367 88483 15170 99456 18185 11427 21583 22372 29980 33304 42964
!44253 51964 55219 59158 66207 60960 77185 62963 88182 60965 99174 58373 11169
!59179 22167 59181 33166 589// 44164 579// 55160 591// 66155 567// 77150 565//
!88148 557// 99139 569// 11135 559// 22116 575// 33113 567// 44106 563// 55103
!577// 66102 567// 31313 58208 82326 41414 21701
! PPBB 56008 72694 90012 29006 27508 26008 90346 26007 28006 24509 90789 21517
!20519 22522 91246 25525 25027 25030 92057 24537 24540 25046 93017 23543 23045
!23563 939// 22548 9436/ 23554 24034 950// 23040
! TTCC 56001 72694 70864 571// 20017 50078 545// 34004 30407 497// 02010 20675
!471// 05510 10145 339// 07017 88999 77999 31313 58208 82326
!
! If the string 'TTAA' is present, then this coded format is assumed, otherwise
! the 'list' format is assumed.
! The header is read line by line until a valid row of numbers is read with
! columns expected in the order below.
!
!<HTML>
!<TITLE>University of Wyoming - Radiosonde Data</TITLE>
!<LINK REL="StyleSheet" HREF="/resources/select.css" TYPE="text/css">
!<BODY BGCOLOR="white">
!<H2>70273 PANC Anchorage Observations at 00Z 30 May 2017</H2>
!<PRE>
!-----------------------------------------------------------------------------
! PRES HGHT TEMP DWPT RELH MIXR DRCT SKNT THTA THTE THTV
! hPa m C C % g/kg deg knot K K K
!-----------------------------------------------------------------------------
! 1008.0 40 7.6 4.9 83 5.41 240 3 280.1 295.2 281.0
! 1003.0 91 5.8 3.8 87 5.03 243 4 278.7 292.7 279.6
! 1000.0 121 5.8 3.7 86 5.01 245 4 278.9 292.9 279.8
!:
!:
!:
! 13.3 29566 -38.5 -66.3 4 0.38 65 16 805.6 809.9 805.8
! 12.2 30175 -37.1 -66.2 3 0.42 110 23 831.2 836.1 831.4
! 11.7 30480 -36.4 -66.2 3 0.45 115 12 844.4 849.5 844.6
! 11.2 30785 -35.7 -66.1 3 0.47 90 20 857.7 863.2 857.9
! 10.7 31072 -35.1 -66.1 3 0.49 870.4 876.3 870.6
!</PRE><H3>Station information and sounding indices</H3><PRE>
! Station identifier: PANC
! Station number: 70273
! Observation time: 170530/0000
! Station latitude: 61.16
! Station longitude: -150.01
! Station elevation: 40.0
! Showalter index: 1.68
! Lifted index: 2.33
! LIFT computed using virtual temperature: 2.32
! SWEAT index: 228.00
! K index: 21.60
! Cross totals index: 28.10
! Vertical totals index: 29.80
! Totals totals index: 57.90
! Convective Available Potential Energy: 0.00
! CAPE using virtual temperature: 0.00
! Convective Inhibition: 0.00
! CINS using virtual temperature: 0.00
! Bulk Richardson Number: 0.00
! Bulk Richardson Number using CAPV: 0.00
! Temp [K] of the Lifted Condensation Level: 274.14
!Pres [hPa] of the Lifted Condensation Level: 934.80
! Mean mixed layer potential temperature: 279.49
! Mean mixed layer mixing ratio: 4.44
! 1000 hPa to 500 hPa thickness: 5289.00
!Precipitable water [mm] for entire sounding: 11.68
!</PRE>
write(linebuffer080 ,'(a80)') repeat(' ',80)
write(linebuffer080_2,'(a80)') repeat(' ',80)
write(linebuffer080_3,'(a80)') repeat(' ',80)
write(linebuffer080_4,'(a80)') repeat(' ',80)
MR_SndVarsName(:) = " "
MR_SndVarsName( 0) = "P "
MR_SndVarsName( 1) = "H "
MR_SndVarsName( 2) = "U "
MR_SndVarsName( 3) = "V "
MR_SndVarsName( 4) = "W "
MR_SndVarsName( 5) = "T "
! Note that there is a deviation here from the master list: Met_var_NC_names
! Met_var:: 6 is for 3d pressure
! Met_var:: 7 is for pressure vert. velocity
MR_SndVarsName( 6) = "Wsp"
MR_SndVarsName( 7) = "Wdr"
MR_SndVarsName(30) = "RH "
MR_SndVarsName(31) = "SH "
MR_SndVarsName(32) = "QL "
MR_SndVarsName(33) = "QI "
MR_SndVarsName(44) = "P0 "
IsCustVarOrder = .false.
IsWindDirectSpeed = .true.
Met_dim_IsAvailable = .false.
Met_var_IsAvailable = .false.
MR_nstat = min(MR_nSnd_Locs,MR_nstat)
! HFS: Break up reader into a diagnostic (type of ASCII file), then subroutines
! for the individual types (ASCII column, Radiosonde: list, raw, RUC)
! This is the start of a huge if statement that determines the type of ASCII input
! and loads the data. All data for all time steps are loaded here into the
! variable MR_SndVars_metP(MR_nSnd_Locs,MR_Snd_nt_fullmet,MR_Snd_nvars,nrows)
! The if statement sorts into the following cases:
! 1 : ASCII sonde file with either 3, 5, or a custom number of columns of data
! 2 : radio sonde data either the clear text variety or the WMO/GTS encoded
if(MR_iwind.eq.1.and.MR_iwindformat.eq.1)then
! We are reading just one windfile with the following format
! If three values per line:
! Elevation(m) , Velocity(m/s) , Direction(degree E of N)
! If five values per line:
! Elevation(m) , Velocity(m/s) , Direction(degree E of N) , Pressure(hPa) ,
! Temperature(C)
! Or a custom list of variables if specified on line 2
! x and y fullmet array will just be the coordinates in the order listed
allocate(x_fullmet_sp(MR_nSnd_Locs))
allocate(y_fullmet_sp(MR_nSnd_Locs))
allocate(MR_dx_met(MR_nSnd_Locs))
allocate(MR_dy_met(MR_nSnd_Locs))
allocate(Snd_idx(MR_nSnd_Locs))
! The number of time steps/file is fixed to 1
nt_fullmet = 1
! There may be more windfiles than time steps if there are multiple sonde locations,
! but we will need a time stamp on each windfile
allocate(MR_windfile_starthour(MR_iwindfiles))
allocate(MR_windfile_stephour(MR_iwindfiles,nt_fullmet))
MR_windfile_starthour = 0.0_dp ! This is the initialization; will be set below
MR_windfile_stephour = 0.0_dp ! This will be the final value since all files
! have one step and so no offset.
MR_windfiles_nt_fullmet(:) = nt_fullmet
do itime = 1,MR_Snd_nt_fullmet
do iloc = 1,MR_nSnd_Locs
iw_idx = (itime-1)*MR_nSnd_Locs + iloc
do io=1,MR_nio;if(VB(io).le.verbosity_info)then
write(outlog(io),*)"Opening sonde file ",iw_idx,&
trim(adjustl(MR_windfiles(iw_idx)))
endif;enddo
open(unit=fid,file=trim(adjustl(MR_windfiles(iw_idx))), status='old',action='read',err=1971)
! skip over first line (Comment line, maybe the location name)
read(fid,*,iostat=iostatus,iomsg=iomessage)linebuffer080
if(iostatus.ne.0) call MR_FileIO_Error_Handler(iostatus,linebuffer080,iomessage)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Reading L2: time(hr) nlev [ncol] [ivar(ncol)]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
read(fid,'(a80)',iostat=iostatus,iomsg=iomessage)linebuffer080
if(iostatus.ne.0) call MR_FileIO_Error_Handler(iostatus,linebuffer080,iomessage)
! Assume we can read at least two values (a real and an integer)
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage) rvalue1, ivalue1
if(iostatus.ne.0)then
do io=1,MR_nio;if(VB(io).le.verbosity_error)then
if(iostatus.lt.0)then
write(errlog(io),*)'MR ERROR: EOR encountered'
else
write(errlog(io),*)'MR ERROR: Error reading line 2 of ASCII windfile'
write(errlog(io),*)' Expecting to read: rvalue1, ivalue1'
write(errlog(io),*)' From the following line from the file: '
write(errlog(io),*)linebuffer080
write(errlog(io),*)'MR System Message: ',trim(adjustl(iomessage))
endif
endif;enddo
stop 1
endif
WindTime = rvalue1
MR_windfile_starthour(iw_idx) = WindTime
nlev = ivalue1
! Try for three values [ncol]
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage) rvalue1, ivalue1, ivalue2
if(iostatus.eq.0)then
! Success: third value is the number of variables
! Note: We need at least 5 variables for P,H,U,V,T
! First check if this is the first file read, otherwise do not allocate
if(iw_idx.eq.1)then
IsCustVarOrder = .true.
MR_Snd_nvars = max(5,ivalue2)
allocate(MR_SndVarsID(MR_Snd_nvars)) ! This is the storage oder
MR_SndVarsID(1) = 0 ! P
MR_SndVarsID(2) = 1 ! H
MR_SndVarsID(3) = 2 ! U
MR_SndVarsID(4) = 3 ! V
MR_SndVarsID(5) = 5 ! T
! There might be more columns to map
allocate(SndColReadOrder(MR_Snd_nvars)) ! This is the read oder
endif
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage) &
rvalue1,ivalue1, ivalue2, SndColReadOrder(1:MR_Snd_nvars)
do io=1,MR_nio;if(VB(io).le.verbosity_info)then
write(outlog(io),*)" Parsing sonde info line as:"
write(outlog(io),*)" wind time: ",rvalue1
write(outlog(io),*)" nlevels : ",ivalue1
write(outlog(io),*)" nvars : ",ivalue2
write(outlog(io),*)" ivars(1:nvar) : ",SndColReadOrder(1:MR_Snd_nvars)
write(outlog(io),*)"1-d ASCII file should contain ",ivalue2," columns"
endif;enddo
do iv = 1,MR_Snd_nvars
do io=1,MR_nio;if(VB(io).le.verbosity_info)then
write(outlog(io),*)" Column ",iv,MR_SndVarsName(SndColReadOrder(iv))
endif;enddo
if (SndColReadOrder(iv).eq.0)then
Met_dim_IsAvailable(2) = .true. ! P
elseif (SndColReadOrder(iv).eq.1)then
Met_var_IsAvailable(1) = .true. ! GPH
elseif (SndColReadOrder(iv).eq.2)then
Met_var_IsAvailable(2) = .true. ! U
elseif (SndColReadOrder(iv).eq.3)then
Met_var_IsAvailable(3) = .true. ! V
elseif (SndColReadOrder(iv).eq.4)then
Met_var_IsAvailable(4) = .true. ! W
Met_var_IsAvailable(7) = .true. ! VVP
elseif (SndColReadOrder(iv).eq.5)then
Met_var_IsAvailable(5) = .true. ! T
elseif (SndColReadOrder(iv).eq.6)then
Met_var_IsAvailable(8) = .true. ! Wsp
elseif (SndColReadOrder(iv).eq.7)then
Met_var_IsAvailable(9) = .true. ! Wdr
endif
enddo
! Now checking how velocities are provided. If coordinates are given, they take
! precedence and are used, with direction and speed ignored. Otherwise, use
! direction and speed.
if(Met_var_IsAvailable(2).and.Met_var_IsAvailable(3))then
IsWindDirectSpeed = .false.
if(Met_var_IsAvailable(2).and.Met_var_IsAvailable(8))then
do io=1,MR_nio;if(VB(io).le.verbosity_info)then
write(outlog(io),*)"MR WARNING: Both U/V and Speed/Direction provided."
write(outlog(io),*)" Ignoring Speed/Direction"
endif;enddo
endif
else
IsWindDirectSpeed = .true.
endif
if (.not.(Met_dim_IsAvailable(2).or.Met_var_IsAvailable(1)))then
do io=1,MR_nio;if(VB(io).le.verbosity_error)then
write(errlog(io),*)"MR ERROR: No height variable given"
endif;enddo
stop 1
elseif(.not.(Met_var_IsAvailable(2).and.Met_var_IsAvailable(3)).and. &
.not.(Met_var_IsAvailable(8).and.Met_var_IsAvailable(9)))then
do io=1,MR_nio;if(VB(io).le.verbosity_error)then
write(errlog(io),*)"MR ERROR: No U/V or Wsp/Wdr variable given"
endif;enddo
stop 1
endif
else ! end of branch where 3 values are read in
! If no list of variables is provided, we will still need a list up to 5
if(iw_idx.eq.1)then
MR_Snd_nvars = 5
allocate(MR_SndVarsID(MR_Snd_nvars))
! This maps the column index of MR_SndVars_metP to ivar
MR_SndVarsID(1) = 0 ! P
MR_SndVarsID(2) = 1 ! H
MR_SndVarsID(3) = 2 ! U
MR_SndVarsID(4) = 3 ! V
MR_SndVarsID(5) = 5 ! T
endif
endif
! Now we know what the columns will mean (3-col, 5-col, custom)
! This variable will hold all the sonde data, up to MAX_ROWS rows
! The order of the columns will be P, H, U, V, T + extras if given
if(iw_idx.eq.1)then
allocate(MR_SndVars_metP(MR_nSnd_Locs,MR_Snd_nt_fullmet,MR_Snd_nvars,MAX_ROWS))
allocate(MR_Snd_np_fullmet(MR_nSnd_Locs,MR_Snd_nt_fullmet))
MR_SndVars_metP = 0.0_sp
MR_Snd_np_fullmet = 0
endif
MR_Snd_np_fullmet(iloc,itime) = nlev
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Reading L3: x/lon y/lat [LLflag] [proj flags]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
read(fid,'(a80)',iostat=iostatus,iomsg=iomessage)linebuffer080
if(iostatus.ne.0) call MR_FileIO_Error_Handler(iostatus,linebuffer080,iomessage)
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)&
rvalue1, rvalue2 ! These are the coordinates of the sonde point
if(iostatus.ne.0)then
do io=1,MR_nio;if(VB(io).le.verbosity_error)then
if(iostatus.lt.0)then
write(errlog(io),*)'MR ERROR: EOF encountered'
else
write(errlog(io),*)'MR ERROR: Error reading line 3 of ASCII windfile'
write(errlog(io),*)' Expecting to read: rvalue1, rvalue2'
write(errlog(io),*)' From the follwing line from the file: '
write(errlog(io),*)linebuffer080
endif
endif;enddo
stop 1
endif
x_fullmet_sp(iloc) = rvalue1
y_fullmet_sp(iloc) = rvalue2
! Try for three values
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage) rvalue1,rvalue2, ivalue1
if(iostatus.eq.0)then
! A third parameter is present: first value of projection line
Snd_Have_Coord = .true.
ilatlonflag = ivalue1
else
! A third parameter is not present: assume met coordinates are Lon/Lat
Snd_Have_Coord = .false.
ilatlonflag = 1
IsLatLon_MetGrid = .true.
endif
if(Snd_Have_Coord)then
if(ilatlonflag.eq.1)then
! We are using a Lon/Lat grid. No need to read anything
! else on this line.
IsLatLon_MetGrid = .true.
IsGlobal_MetGrid = .false.
IsRegular_MetGrid = .false.
! The rest of this is ignored
Met_iprojflag = 4
Met_lam0 = 265.0_8
Met_phi0 = 25.0_8
Met_phi1 = 25.0_8
Met_phi2 = 25.0_8
Met_k0 = 0.933_8
Met_Re = 6371.229_8
else
! Grid in not Lon/Lat
! First try to read the projection code
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage) rvalue1,rvalue2, &
ivalue1, ivalue2
if(iostatus.eq.0)then
Met_iprojflag = ivalue2
if(Met_iprojflag.ne.0)then
! we have a geographic projection
! Try to read the full projection line
! We know the projection flag will be 0 since this is the branch for
! that so we know ' 0 ' will be in the string, but the first or
! second coordinate could have this string too. Need a better way.
! ToDo: Need to compress repeated white spaces then look for the second space
indx1 = index(linebuffer080,' 0 ')
indx2 = index(linebuffer080,'#') ! Check for the comment marker
indx3 = len_trim(linebuffer080) ! get length of string
if(indx2.gt.0)then
linebuffer080_2(1:indx1+indx2) = linebuffer080(indx1+1:indx2-1)
else
linebuffer080_2(1:indx1+1+indx3) = linebuffer080(indx1+1:indx3)
endif
call PJ_Set_Proj_Params(linebuffer080_2)
IsLatLon_MetGrid = .false.
IsGlobal_MetGrid = .false.
IsRegular_MetGrid = .false.
Met_iprojflag = PJ_iprojflag
Met_k0 = PJ_k0
Met_Re = PJ_Re
Met_lam0 = PJ_lam0
Met_lam1 = PJ_lam1
Met_lam2 = PJ_lam2
Met_phi0 = PJ_phi0
Met_phi1 = PJ_phi1
Met_phi2 = PJ_phi2
endif
else
! If we have a cartesian grid, but no projection code,
! assume no geographic projection
Met_iprojflag=0
endif
endif
else
! No projection info given, assume Lon/Lat
IsLatLon_MetGrid = .true.
endif
! Finished projection parameters,
! Allocating temporary space for wind data
allocate( WindVelocity(nlev)); WindVelocity = 0.0_sp
allocate(WindDirection(nlev)); WindDirection = 0.0_sp
!Read elevation (m), speed, direction at each level
! Speed is given in m/s (multiply by 0.514444444 to convert
! from knots to m/s)
! Direction is given as degrees east of north and specifies the
! direction FROM which the wind is blowing.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Reading L4-EOF: ncols of data in nlev rows
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do il=1,nlev
read(fid,'(a80)',iostat=iostatus,iomsg=iomessage)linebuffer080
if(iostatus.ne.0) call MR_FileIO_Error_Handler(iostatus,linebuffer080,iomessage)
! For the first line, we need to determine the number of columns.
if (il.eq.1)then
! Verify we can at least read 3
ncols = 0
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)rvalues(1:3)
if(iostatus.eq.0)then
ncols=3
! Now try 4
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)rvalues(1:4)
if(iostatus.eq.0)then
ncols=4
! Now try 5
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)rvalues(1:5)
if(iostatus.eq.0)then
ncols=5
! Now try 6
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)rvalues(1:6)
if(iostatus.eq.0)then
ncols=6
! Now try 7
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)rvalues(1:7)
if(iostatus.eq.0)then
ncols=7
! Now try 8
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)rvalues(1:8)
if(iostatus.eq.0)then
ncols=8
! Now try 9
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)rvalues(1:9)
if(iostatus.eq.0)then
ncols=9
! Now try 10
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage)rvalues(1:10)
if(iostatus.eq.0)then
ncols=10
endif
endif
endif
endif
endif
endif
endif
else
do io=1,MR_nio;if(VB(io).le.verbosity_error)then
write(errlog(io),*)&
"MR ERROR: ASCII wind files must have at least 3 column of data."
endif;enddo
stop 1
endif
endif ! il.eq.1
rvalues(:) = -1.99_sp
! read ncols of data on this row
read(linebuffer080,*,iostat=iostatus,iomsg=iomessage) rvalues(1:ncols)
if(.not.IsCustVarOrder.and.ncols.eq.3)then
! If this is a 3-column file without the custom ordering, read as
! Altitude (m), wind speed (m/s), wind direction (deg E of N)
IsWindDirectSpeed = .true.
MR_SndVars_metP(iloc,itime,2,il) = rvalues(1)*1.0e-3_sp ! convert to km
WindVelocity(il) = rvalues(2)
WindDirection(il) = rvalues(3)
! Calculate pressure in Pa and T in K
MR_SndVars_metP(iloc,itime,1,il) = MR_Pres_US_StdAtm(MR_SndVars_metP(iloc,itime,2,il))&
* 100.0_sp
MR_SndVars_metP(iloc,itime,5,il) = MR_Temp_US_StdAtm(MR_SndVars_metP(iloc,itime,2,il))
elseif(.not.IsCustVarOrder.and.ncols.eq.5)then
! If this is a 5-column file without the custom ordering, read as
! Recall the first five columns are P,H,U,V,T
IsWindDirectSpeed = .true.
Snd_Have_PT = .true.
! Five values were successfully read, interpret as:
! Altitude (m), wind speed (m/s), wind direction (deg E of N)
! Pressure (hPa), Temperature (C)
MR_SndVars_metP(iloc,itime,2,il) = rvalues(1)*1.0e-3_sp ! convert to km
WindVelocity(il) = rvalues(2)
WindDirection(il) = rvalues(3)
if(iw_idx.eq.1.and.& ! For the first file
il.eq.1.and. & ! and the first level
rvalues(4).gt.1500.0_sp)then ! test pressure value
! Check if pressure is greater than the expected 1013 hPa.
! If so, assume pressure is in Pa
In_hPa = .false.
endif
if(In_hPa)then
MR_SndVars_metP(iloc,itime,1,il) = rvalues(4) * 100.0_sp ! convert to Pa
else
MR_SndVars_metP(iloc,itime,1,il) = rvalues(4) ! already in Pa
endif
MR_SndVars_metP(iloc,itime,5,il) = rvalues(5) + 273.0_sp ! convert to K
elseif(IsCustVarOrder)then
! This is the customized list of data columns
do ic=1,ncols
! Unfortunately, we currently need to map the variable ID to the index
if(SndColReadOrder(ic).eq.0)then ! pressure
if(iw_idx.eq.1.and.& ! For the first file
il.eq.1.and. & ! and the first level
rvalues(ic).gt.1500.0_sp)then ! test pressure value
In_hPa = .false.
else
In_hPa = .true.
endif
if(In_hPa)then
MR_SndVars_metP(iloc,itime,1,il) = rvalues(ic) * 100.0_sp
else
MR_SndVars_metP(iloc,itime,1,il) = rvalues(ic)
endif
elseif(SndColReadOrder(ic).eq.1)then ! Altitude (m)
MR_SndVars_metP(iloc,itime,2,il) = rvalues(ic)*1.0e-3_sp ! convert to km
elseif(SndColReadOrder(ic).eq.2)then ! U (m/s)
MR_SndVars_metP(iloc,itime,3,il) = rvalues(ic)
elseif(SndColReadOrder(ic).eq.3)then ! V (m/s)
MR_SndVars_metP(iloc,itime,4,il) = rvalues(ic)
!elseif(SndColReadOrder(ic).eq.4)then ! Vert velocity (m/s)
! MR_SndVars_metP(iloc,itime,?,il) = rvalues(i)
elseif(SndColReadOrder(ic).eq.5)then ! Temperature (C)
MR_SndVars_metP(iloc,itime,5,il) = rvalues(ic) + 273.0_sp ! convert to K
elseif(SndColReadOrder(ic).eq.6)then ! Wind speed (m/s)
WindVelocity(il) = rvalues(ic)
elseif(SndColReadOrder(ic).eq.7)then ! Wind direction (from deg E of N)
WindDirection(il) = rvalues(ic)
else
do io=1,MR_nio;if(VB(io).le.verbosity_info)then
write(errlog(io),*)&
"MR WARNING: Ignoring data for variable ",MR_SndVarsName(SndColReadOrder(ic))
endif;enddo
endif
enddo
endif
! Finished interpreting the row of data
if(IsWindDirectSpeed)then
! If the wind data is provided as speed and direction, break it into U,V components
MR_SndVars_metP(iloc,itime,3,il) = &
real(WindVelocity(il)*sin(pi + DEG2RAD*WindDirection(il)),kind=sp)
MR_SndVars_metP(iloc,itime,4,il) = &
real(WindVelocity(il)*cos(pi + DEG2RAD*WindDirection(il)),kind=sp)
else
! Otherwise, calculate direction and speed from the components
WindVelocity(il) = &
MR_SndVars_metP(iloc,itime,3,il)*MR_SndVars_metP(iloc,itime,3,il) + &
MR_SndVars_metP(iloc,itime,4,il)*MR_SndVars_metP(iloc,itime,4,il)
WindVelocity(il) = sqrt(WindVelocity(il))
WindDirection(il) = &
-real(atan2(MR_SndVars_metP(iloc,itime,4,il),MR_SndVars_metP(iloc,itime,3,il))/DEG2RAD,kind=sp) + &
90.0_sp
endif
!Met_dim_IsAvailable(1) = .true. ! Time
Met_dim_IsAvailable(2) = .true. ! P
!Met_dim_IsAvailable(3) = .true. ! y
!Met_dim_IsAvailable(4) = .true. ! x
Met_var_IsAvailable(1) = .true. ! GPH
Met_var_IsAvailable(2) = .true. ! U
Met_var_IsAvailable(3) = .true. ! V
Met_var_IsAvailable(5) = .true. ! T
if(il.eq.1)then
do io=1,MR_nio;if(VB(io).le.verbosity_info)then
write(outlog(io),*)"Sonde values loaded."
write(outlog(io),203)MR_SndVarsName(1:7)
endif;enddo
endif
do io=1,MR_nio;if(VB(io).le.verbosity_info)then
write(outlog(io),204)il,real(MR_SndVars_metP(iloc,itime,1:5,il),kind=4),&
WindVelocity(il),WindDirection(il)
endif;enddo
203 format(3x,'Level',6x,a3,'(Pa)', &
7x,a3,'(km)', &
3x,a3,'(m/s)', &
2x,a3,'(m/s)', &
1x,a3,'(K)', &
3x,a3,'(m/s)', &
2x,a3,'(degree E of N)')
204 format(6x,i2,5x,7f10.3)
enddo ! il=1,nlev
close(fid)
! Finished reading all the data for this file
deallocate( WindVelocity)
deallocate(WindDirection)
enddo ! iloc = 1,MR_nSnd_Locs
enddo ! itime = 1,MR_Snd_nt_fullmet
! Finished reading all the data of all the files
! Here we look for the highest pressure value (lowest altitude) of all the pressure
! values to be used for setting the master pressure array (p_fullmet_sp)
p_maxtop = 0.0_sp
p_tidx = 0
p_lidx = 0
do itime = 1,MR_Snd_nt_fullmet
do iloc = 1,MR_nSnd_Locs
p_top = MR_SndVars_metP(iloc,itime,1,MR_Snd_np_fullmet(iloc,itime))
if(p_top.gt.p_maxtop)then
p_maxtop = MR_Snd_np_fullmet(iloc,itime)
p_tidx = itime
p_lidx = iloc
endif
enddo
enddo
np_fullmet = MR_Snd_np_fullmet(p_lidx,p_tidx)
!np_fullmet_Vz = np_fullmet
!np_fullmet_RH = np_fullmet
allocate(p_fullmet_sp(np_fullmet))
!allocate(p_fullmet_Vz_sp(np_fullmet_Vz))
!allocate(p_fullmet_RH_sp(np_fullmet_RH))
p_fullmet_sp(1:np_fullmet)=MR_SndVars_metP(p_lidx,p_tidx,1,1:np_fullmet)
!p_fullmet_Vz_sp = p_fullmet_sp
!p_fullmet_RH_sp = p_fullmet_sp
MR_Max_geoH_metP_predicted = MR_Z_US_StdAtm(p_fullmet_sp(np_fullmet)/100.0_sp)
allocate(z_approx(np_fullmet))
do k=1,np_fullmet
! Calculate heights for US Std Atmos while pressures are still in mbars
! or hPa
z_approx(k) = MR_Z_US_StdAtm(p_fullmet_sp(k))
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
elseif(MR_iwind.eq.1.and.MR_iwindformat.eq.2)then
! We are reading radiosonde data from http://weather.uwyo.edu/ or https://ruc.noaa.gov/raobs/
! First, these are always given in lon/lat, earth-relative coordinates:
IsLatLon_MetGrid = .true.
IsGlobal_MetGrid = .false.
IsRegular_MetGrid = .false.
! The rest of this is ignored
Met_iprojflag = 4
Met_lam0 = 265.0_8
Met_phi0 = 25.0_8
Met_phi1 = 25.0_8
Met_phi2 = 25.0_8
Met_k0 = 0.933_8
Met_Re = 6371.229_8
! Load the internal database of radiosonde stations
call MR_Load_Radiosonde_Station_Data
! Allocate arrays for iGridCode points at iwindfiles times, and MAX_ROWS levels
! x and y fullmet array will just be the coordinates in the order listed
allocate(x_fullmet_sp(MR_nSnd_Locs))
allocate(y_fullmet_sp(MR_nSnd_Locs))
allocate(MR_dx_met(MR_nSnd_Locs))
allocate(MR_dy_met(MR_nSnd_Locs))
allocate(Snd_idx(MR_nSnd_Locs))
! The number of time steps/file is fixed to 1
nt_fullmet = 1
! There may be more windfiles than time steps if there are multiple sonde locations,
! but we will need a time stamp on each windfile
allocate(MR_windfile_starthour(MR_iwindfiles))
allocate(MR_windfile_stephour(MR_iwindfiles,nt_fullmet))
MR_windfile_starthour = 0.0_dp ! This is the initialization; will be set below
MR_windfile_stephour = 0.0_dp ! This will be the final value since all files
! have one step and so no offset.
MR_windfiles_nt_fullmet(:) = nt_fullmet
do itime = 1,MR_Snd_nt_fullmet
do iloc = 1,MR_nSnd_Locs
iw_idx = (itime-1)*MR_nSnd_Locs + iloc
do io=1,MR_nio;if(VB(io).le.verbosity_info)then
write(outlog(io),*)"Opening sonde file ",iw_idx,&
trim(adjustl(MR_windfiles(iw_idx)))
endif;enddo
open(unit=fid,file=trim(adjustl(MR_windfiles(iw_idx))),status='old',action='read',err=1971)
if(iw_idx.eq.1)then
MR_Snd_nvars = 5
allocate(MR_SndVarsID(MR_Snd_nvars)) ! This is the storage order
MR_SndVarsID(1) = 0 ! P
MR_SndVarsID(2) = 1 ! H
MR_SndVarsID(3) = 2 ! U
MR_SndVarsID(4) = 3 ! V
MR_SndVarsID(5) = 5 ! T
! We only allocate 16 rows here because we will only read the TTAA and TTCC mandatory levels
! If you want to use the other information in the radiosonde, shape the data into iwindformat=1
! Note that this might be reduced from 16 if any of the data files do not extend up to 10 mb.
nlev = 16
allocate(MR_SndVars_metP(MR_nSnd_Locs,MR_Snd_nt_fullmet,MR_Snd_nvars,nlev))
allocate(MR_Snd_np_fullmet(MR_nSnd_Locs,MR_Snd_nt_fullmet))
allocate(H_tmp(nlev))
allocate(T_tmp(nlev))
MR_SndVars_metP = 0.0_sp
MR_Snd_np_fullmet(1,1) = 0 ! We initialize this to 0 since we don't know how high
! the sonde went
pres_Snd_tmp(1:nlev) = &
(/1000.0_sp, 925.0_sp, 850.0_sp, 700.0_sp, 500.0_sp, &
400.0_sp, 300.0_sp, 250.0_sp, 200.0_sp, 150.0_sp, &
100.0_sp, 70.0_sp, 50.0_sp, 30.0_sp, 20.0_sp, & ! Note: the TTCC part starts at 70 hPa
10.0_sp/) ! Not all sondes go this high!
endif
! Allocating temporary space for wind data
allocate( WindVelocity(nlev)); WindVelocity = 0.0_sp
allocate(WindDirection(nlev)); WindDirection = 0.0_sp
! First search the whole radiosonde file for the string 'TTAA'. If this
! is present, then assume the file is in FAA604 (WMO/GTS or RAW) format
! TTAA denotes manndatory level data (12 p lev): T,dew-point,gph,w-dir,w-sp
! TTBB denotes significant level data for T,RH (points that help reconstruct profile): we ignore this
! PPBB denotes regional and significant wind levels (also ignored)
! TTCC starts the block that extends TTAA to higher levels if the balloon survives
! TTDD significant level data for pressures less than 100 mb
idx = 0
iostatus = 0
do while (iostatus.ge.0)
read(fid,'(a80)',iostat=iostatus,iomsg=iomessage)linebuffer080
! No aborting on read-error here since we will read this file until the EOF
idx =index(linebuffer080,"TTAA") ! Look for string identifying WMO/GTS or RAW (both have this)
idx2=index(linebuffer080,"TTCC") ! Look for string indicating extended data
if (idx.gt.0) then
IsGTS = .true.
endif
if (idx2.gt.0) then
HasTTCC = .true.
endif
enddo
if(IsGTS)then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Reading data GTS format
! This is either "Text: Raw" from http://weather.uwyo.edu/upperair/sounding.html
!1. TTAA 56001 72694 99009 22868 29006 00137 20667 28507 92801 13662 25507 85505
! or "FAA604 format (WMO/GTS)" from https://ruc.noaa.gov/raobs/
!1. SLEMANSLE
!2. 72694 TTAA 56001 72694 99009 22868 29006 00137 20667 28507
! Note that the format is slightly different so we need to find out
! which it is
GTSstr(1:53)="//////" ! Note: each 'element' of GTSstr(:) is really a 6-char string
rewind(fid)
read(fid,'(a80)',iostat=iostatus,iomsg=iomessage)linebuffer080
if(iostatus.ne.0) call MR_FileIO_Error_Handler(iostatus,linebuffer080,iomessage)
idx =index(linebuffer080,"TTAA")
if(idx.gt.0)then
! 'TTAA' is in the first line, this is a file from Uni. Wyoming
IsRUCNOAA = .false.
! the format for the TTAA block is 3 lines of 13 6-char strings
read(fid,'(a80)',iostat=iostatus,iomsg=iomessage)linebuffer080_2
if(iostatus.ne.0) call MR_FileIO_Error_Handler(iostatus,linebuffer080_2,iomessage)
read(fid,'(a80)',iostat=iostatus,iomsg=iomessage)linebuffer080_3
if(iostatus.ne.0) call MR_FileIO_Error_Handler(iostatus,linebuffer080_3,iomessage)
! the dumstr accommodates the 'TTAA' in the first line with the rest going to GTSstr
read(linebuffer080 ,155,iostat=iostatus,iomsg=iomessage)dumstr1,GTSstr(1:12)
if(iostatus.ne.0)then
do io=1,MR_nio;if(VB(io).le.verbosity_error)then
if(iostatus.lt.0)then
write(errlog(io),*)'MR ERROR: EOR encountered'
else
write(errlog(io),*)'MR ERROR: Error reading GTS format sonde file'
write(errlog(io),*)' Expecting to read: dumstr1,GTSstr(1:12)'
write(errlog(io),*)' with format: 13a6'
write(errlog(io),*)' From the following line from the file: '
write(errlog(io),*)linebuffer080
write(errlog(io),*)'MR System Message: '
write(errlog(io),*)iomessage
endif
endif;enddo
stop 1
endif
! Load the second line to the next 13 elements of GTSstr
read(linebuffer080_2,155,iostat=iostatus,iomsg=iomessage)GTSstr(13:25)
if(iostatus.ne.0)then
do io=1,MR_nio;if(VB(io).le.verbosity_error)then
if(iostatus.lt.0)then
write(errlog(io),*)'MR ERROR: EOR encountered'
else
write(errlog(io),*)'MR ERROR: Error reading GTS format sonde file'
write(errlog(io),*)' Expecting to read: GTSstr(13:25)'
write(errlog(io),*)' with format: 13a6'
write(errlog(io),*)' From the following line from the file: '
write(errlog(io),*)linebuffer080
write(errlog(io),*)'MR System Message: '
write(errlog(io),*)iomessage
endif
endif;enddo
stop 1
endif