forked from hschwaiger-usgs/volcano-ash3d-metreader
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMetReader_ASCII.f90
2672 lines (2572 loc) · 210 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
use projection
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
integer, parameter :: MAX_ROWS = 300 ! maximum number of row of data
integer :: ioerr
integer :: istr1,istr2,istr3
integer :: ic,il,iil,iv
integer :: fid
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
character(len=80) :: linebuffer,linebuffer2,linebuffer3,linebuffer4
character(len=6),dimension(53) :: GTSstr
character(len=6) :: dumstr1,dumstr2,dumstr3,dumstr4
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,SurfTemp,SurfDewPoint,SurfWindDir,SurfWindSpeed
integer :: SurfTemp_int
logical :: In_hPa = .true.
integer :: indx1,indx2
integer :: ncols
logical :: IsWindDirectSpeed
logical :: IsCustVarOrder
integer,dimension(:),allocatable :: SndColReadOrder
character(len=3),dimension(0:50) :: MR_SndVarsName
real(kind=dp) :: HS_hours_since_baseyear ! function that calculates hours
! since base year
integer :: Stat_ID
real(kind=sp) :: Stat_elev
integer :: idx,idx2
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
write(MR_global_production,*)"--------------------------------------------------------------------------------"
write(MR_global_production,*)"---------- MR_Read_Met_DimVars_ASCII_1d ----------"
write(MR_global_production,*)"--------------------------------------------------------------------------------"
!------------------------------------------------------------------------------
! 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>
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 "
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.
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)
! 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))
! 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
write(MR_global_info,*)"Opening sonde file ",iw_idx,&
adjustl(trim(MR_windfiles(iw_idx)))
fid = 127
open(unit=fid,file=trim(adjustl(MR_windfiles(iw_idx))), status='unknown',err=1971)
read(fid,*)!skip over first line
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Reading L2: time(hr) nlev [ncol] [ivar(ncol)]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
read(fid,'(a80)')linebuffer
! Assume we can read at least two values (a real and an integer)
read(linebuffer,*) rvalue1, ivalue1
WindTime = rvalue1
MR_windfile_starthour(iw_idx) = WindTime
nlev = ivalue1
! Try for three values [ncol]
read(linebuffer,*,iostat=ioerr) rvalue1,ivalue1, ivalue2
if(ioerr.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(linebuffer,*,iostat=ioerr) rvalue1,ivalue1, ivalue2, SndColReadOrder(1:MR_Snd_nvars)
write(MR_global_info,*)"1-d ASCII file contains ",ivalue2," columns"
do iv = 1,MR_Snd_nvars
write(MR_global_info,*)" Column ",iv,MR_SndVarsName(SndColReadOrder(iv))
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
elseif (SndColReadOrder(iv).eq.5)then
Met_var_IsAvailable(5) = .true. ! T
elseif (SndColReadOrder(iv).eq.6)then
Met_var_IsAvailable(6) = .true. ! Wsp
elseif (SndColReadOrder(iv).eq.7)then
Met_var_IsAvailable(7) = .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(6).and.Met_var_IsAvailable(7))then
write(MR_global_error,*)"MR WARNING: Both U/V and Speed/Direction provided."
write(MR_global_error,*)" Ignoring Speed/Direction"
endif
else
IsWindDirectSpeed = .true.
endif
if (.not.(Met_dim_IsAvailable(2).or.Met_var_IsAvailable(1)))then
write(MR_global_error,*)"MR ERROR: No height variable given"
stop 1
elseif(.not.(Met_var_IsAvailable(2).and.Met_var_IsAvailable(3)).and. &
.not.(Met_var_IsAvailable(6).and.Met_var_IsAvailable(7)))then
write(MR_global_error,*)"MR ERROR: No U/V or Wsp/Wdr variable given"
stop 1
endif
else
! 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)')linebuffer
read(linebuffer,*,iostat=ioerr) rvalue1,rvalue2 ! These are the coordinates of the sonde point
x_fullmet_sp(iloc) = rvalue1
y_fullmet_sp(iloc) = rvalue2
! Try for three values
read(linebuffer,*,iostat=ioerr) rvalue1,ivalue1, ivalue2
if(ioerr.eq.0)then
! A third parameter is present: first value of projection line
Snd_Have_Coord = .true.
Met_iprojflag = ivalue2
endif
if(Snd_Have_Coord)then
if(Met_iprojflag.eq.1)then
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
! Try to read the projection line
indx1 = index(linebuffer,' 0 ')
indx2 = index(linebuffer,'#')
if(indx2.gt.0)then
call PJ_Set_Proj_Params(linebuffer(indx1:indx2-1))
else
call PJ_Set_Proj_Params(linebuffer(indx1:))
endif
IsLatLon_MetGrid = .false.
IsGlobal_MetGrid = .false.
IsRegular_MetGrid = .false.
Met_iprojflag = PJ_iprojflag
Met_k0 = PJ_k0
Met_Re = PJ_radius_earth
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
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)')linebuffer
! 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(linebuffer,*,iostat=ioerr)rvalues(1:3)
if(ioerr.eq.0)then
ncols=3
! Now try 4
read(linebuffer,*,iostat=ioerr)rvalues(1:4)
if(ioerr.eq.0)then
ncols=4
! Now try 5
read(linebuffer,*,iostat=ioerr)rvalues(1:5)
if(ioerr.eq.0)then
ncols=5
! Now try 6
read(linebuffer,*,iostat=ioerr)rvalues(1:6)
if(ioerr.eq.0)then
ncols=6
! Now try 7
read(linebuffer,*,iostat=ioerr)rvalues(1:7)
if(ioerr.eq.0)then
ncols=7
! Now try 8
read(linebuffer,*,iostat=ioerr)rvalues(1:8)
if(ioerr.eq.0)then
ncols=8
! Now try 9
read(linebuffer,*,iostat=ioerr)rvalues(1:9)
if(ioerr.eq.0)then
ncols=9
! Now try 10
read(linebuffer,*,iostat=ioerr)rvalues(1:10)
if(ioerr.eq.0)then
ncols=10
endif
endif
endif
endif
endif
endif
endif
else
write(MR_global_error,*)&
"MR ERROR: ASCII wind files must have at least 3 column of data."
stop 1
endif
endif ! il.eq.1
rvalues(:) = -1.99_sp
! read ncols of data on this row
read(linebuffer,*,iostat=ioerr) 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
write(MR_global_error,*)&
"MR Warning: Ignoring data for variable ",MR_SndVarsName(SndColReadOrder(ic))
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)
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
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)
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))
write(*,*)k,p_fullmet_sp(k),z_approx(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
! 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))
! 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
write(MR_global_info,*)"Opening sonde file ",iw_idx,&
adjustl(trim(MR_windfiles(iw_idx)))
fid = 127
open(unit=fid,file=trim(adjustl(MR_windfiles(iw_idx))), status='unknown',err=1971)
if(iw_idx.eq.1)then
MR_Snd_nvars = 5
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
! We only allocate 15 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
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
!read(fid,'(a80)',iostat=ioerr)linebuffer
!idx=index(linebuffer,"TTAA")
idx = 0
ioerr = 0
do while (ioerr.ge.0)
read(fid,'(a80)',iostat=ioerr)linebuffer
idx =index(linebuffer,"TTAA")
idx2=index(linebuffer,"TTCC")
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)="//////"
rewind(fid)
read(fid,'(a80)',iostat=ioerr)linebuffer
idx =index(linebuffer,"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=ioerr)linebuffer2
read(fid,'(a80)',iostat=ioerr)linebuffer3
read(linebuffer ,155)dumstr1,GTSstr(1:12) ! the dumstr accommodates the 'TTAA'
read(linebuffer2,155)GTSstr(13:25)
read(linebuffer3,155)GTSstr(26:38)
! Now we need to load the TTCC string blocks if the sonde went high enough
if(HasTTCC)then
idx2=0
do while (idx2.eq.0.and.ioerr.ge.0)
read(fid,'(a80)',iostat=ioerr)linebuffer
idx2=index(linebuffer,"TTCC")
if (idx2.gt.0) then
exit
endif
enddo
read(fid,'(a80)',iostat=ioerr)linebuffer2
read(linebuffer ,155)dumstr1,dumstr2,dumstr3,GTSstr(39:48) ! the dumstr accommodate the first
read(linebuffer2,155)GTSstr(49:53) ! three blocks
endif
else
! TTAA not in the first line, try the second
read(fid,'(a80)',iostat=ioerr)linebuffer
idx =index(linebuffer,"TTAA")
if(idx.gt.0)then
! this is a file from NOAA Rapid Update Cycle ruc.noaa.gov
IsRUCNOAA = .true.
! the format for the TTAA block is 4 lines of 10 6-char strings
read(fid,'(a80)',iostat=ioerr)linebuffer2
read(fid,'(a80)',iostat=ioerr)linebuffer3
read(fid,'(a80)',iostat=ioerr)linebuffer4
read(linebuffer ,154)dumstr1,dumstr2,GTSstr(1:8) ! here we need the space
! for the repeat StatID and TTAA
read(linebuffer2,154)GTSstr( 9:18)
read(linebuffer3,154)GTSstr(19:28)
read(linebuffer4,154)GTSstr(29:38)
! Now we need to load the TTCC string blocks if the sonde went high enough
if(HasTTCC)then
idx2=0
do while (idx2.eq.0.and.ioerr.ge.0)
read(fid,'(a80)',iostat=ioerr)linebuffer
idx2=index(linebuffer,"TTCC")
if (idx2.gt.0) then
exit
endif
enddo
read(fid,'(a80)',iostat=ioerr)linebuffer2
read(linebuffer ,154)dumstr1,dumstr2,dumstr3,dumstr4,GTSstr(39:44)
read(linebuffer2,154)GTSstr(45:53)
endif
else
! This file has a 'TTAA' in it, but is not one of the expected
! formats (first or second line). Abort.
write(MR_global_error,*)"MR ERROR: Cannot determine the format of the"
write(MR_global_error,*)" radiosonde data."
stop 1
endif
endif
! Now interpret the strings
! Block B: day/time
read(GTSstr(1)(1:2),*)DayOfMonth
DayOfMonth = DayOfMonth-50
read(GTSstr(1)(3:4),*)SndHour
MR_windfile_starthour(iw_idx) = HS_hours_since_baseyear(2018,6,DayOfMonth,&
real(SndHour,kind=dp),MR_BaseYear,MR_useLeap)
! Block A: station identifier
read(GTSstr(2),*)Stat_ID
call MR_Get_Radiosonde_Station_Coord(Stat_ID, &
y_fullmet_sp(iloc),x_fullmet_sp(iloc),&
Stat_elev)
! Block C: surface pressure
if(GTSstr(3)(1:1).ne.'/')then
read(GTSstr(3)(1:2),*)dum_int ! should be 99 indicating surface
else
dum_int = -9999
endif
if(GTSstr(3)(3:3).ne.'/')then
read(GTSstr(3)(3:5),*)dum_int ! surface pressure in mb
else
dum_int = -9999
endif
SurfPres = real(dum_int,kind=sp)
! Block D: temperature/dew-point for last pressure
if(GTSstr(4)(1:1).ne.'/')then
read(GTSstr(4)(1:3),*)SurfTemp_int ! temperature is characters 1-3, dew-point 4-5
else
SurfTemp_int = -9999
endif
if (mod(SurfTemp_int,2).eq.0)then
SurfTemp = 0.1_sp*real(SurfTemp_int,kind=sp)
else
SurfTemp = -0.1_sp*real(SurfTemp_int,kind=sp)
endif
SurfTemp = SurfTemp
if(GTSstr(4)(4:4).ne.'/')then
read(GTSstr(4)(4:5),*)dum_int
else
dum_int = -9999
endif
SurfDewPoint = real(dum_int,kind=sp)
! Block E: Wind direction and speed for last pressure
if(GTSstr(5)(1:1).ne.'/')then
read(GTSstr(5)(1:3),*)dum_int ! Wind direction is char 1-3
else
dum_int = -9999
endif
SurfWindDir = real(dum_int,kind=sp)
if(GTSstr(5)(4:4).ne.'/')then
read(GTSstr(5)(4:5),*)dum_int ! Wind speed (knts) is char 4-5
else
dum_int = -9999
endif
SurfWindSpeed = real(dum_int,kind=sp)
if(HasTTCC)then
ulev = nlev
else
ulev = 11
endif
MR_Snd_np_fullmet(iloc,itime) = ulev
do il = 1,ulev
istr1 = 6 + (il-1)*3 ! block for pressure/height
istr2 = 6 + (il-1)*3 +1 ! block for temperature/dew-point
istr3 = 6 + (il-1)*3 +2 ! block for wind direction/speed
! Block 1 of the 3-block level data (pressure, height)
if(GTSstr(istr1)(1:1).ne.'/')then
read(GTSstr(istr1)(1:2),*)dum_int ! this is a pressure indicator that we could double-check
else
dum_int = -9999
endif
MR_SndVars_metP(iloc,itime,1,il) = pres_Snd_tmp(il)
if(GTSstr(istr1)(3:3).ne.'/')then
read(GTSstr(istr1)(3:5),*)H_tmp(il) ! height a.s.l
else
H_tmp(il) = -9999
endif
! Block 2 of the 3-block level data (temperature, dew point (not read))
if(GTSstr(istr2)(1:1).ne.'/')then
read(GTSstr(istr2)(1:3),*)T_tmp(il) ! temperature in C is characters 1-3, dew-point 4-5
else
T_tmp(il) = -9999
endif
! Block 3 of the 3-block level data (wind direction, wind speed)
if(GTSstr(istr3)(1:1).ne.'/')then
read(GTSstr(istr3)(1:3),*)dum_int ! Wind direction is char 1-3
else
dum_int = -9999
endif
WindDirection(il) = real(dum_int,kind=sp)
if(GTSstr(istr3)(4:4).ne.'/')then
read(GTSstr(istr3)(4:5),*)dum_int ! Wind speed (knts) is char 4-5
else
dum_int = -9999
endif
WindVelocity(il) = real(dum_int,kind=sp)
enddo
!--------------------------------
! Now convert to the proper units.
! Pressure: No error-checking
MR_SndVars_metP(iloc,itime,1,:) = MR_SndVars_metP(iloc,itime,1,:) * 100.0_sp ! convert to Pa
! Height: also assume valid
! do the block-specific adjustments to get the correct m blocks for 1000 mb
if(H_tmp(1).ge.500)then ! kludge for negative heights
MR_SndVars_metP(iloc,itime,2,1) = -1.0_sp * real(H_tmp(1)-500,kind=sp)
else
MR_SndVars_metP(iloc,itime,2,1) = real(H_tmp(1),kind=sp)
endif
! The heights are adjusted with a multiplicative and additive factor. If the
! three-digit value is less then the previous, then we are in the next scaling
! bracket
scl_idx = 1
scl_m(1:7) = (/1.0_sp, 1.0_sp, 1.0_sp, 10.0_sp, &
10.0_sp, 10.0_sp, 10.0_sp/)
scl_a(1:7) = (/0.0_sp, 1000.0_sp, 3000.0_sp, 0.0_sp, &
10000.0_sp, 20000.0_sp, 30000.0_sp/)
do il = 2,ulev
if (H_tmp(il).lt.H_tmp(il-1))then
! check if we need to increment the scaling bracket
scl_idx = scl_idx + 1
endif
if (il.eq.5)scl_idx=4 ! for the 500 mb level, make sure we are using decimeters
MR_SndVars_metP(iloc,itime,2,il) = scl_a(scl_idx) + scl_m(scl_idx)* &
real(H_tmp(il),kind=sp)
enddo
MR_SndVars_metP(iloc,itime,2,:) = MR_SndVars_metP(iloc,itime,2,:)*1.0e-3_sp ! convert m to km
! Temperature: do error checking
do il = 1,ulev
if(T_tmp(il).ne.-9999)then
if (mod(dum_int,2).eq.0)then
MR_SndVars_metP(iloc,itime,5,il) = 0.1_sp*real(T_tmp(il),kind=sp)
else
MR_SndVars_metP(iloc,itime,5,il) = -0.1_sp*real(T_tmp(il),kind=sp)
endif
else
MR_SndVars_metP(iloc,itime,5,il) = real(T_tmp(il),kind=sp)
endif
enddo
do il = 1,ulev
! Set NaN's below the surface to the surface value
if(pres_Snd_tmp(il).gt.SurfPres.or.il.eq.1)then
if(T_tmp(il).eq.-9999) &
MR_SndVars_metP(iloc,itime,5,il) = SurfTemp
else
! For NaN's above the surface, interpolate to the next valid value
if(T_tmp(il).eq.-9999)then
il_bot=il-1
do ill=il,ulev
if(T_tmp(ill).ne.-9999)then
il_top = ill
exit
endif
enddo
zhere=MR_SndVars_metP(iloc,itime,2,il)
zbot=MR_SndVars_metP(iloc,itime,2,il_bot)
ztop=MR_SndVars_metP(iloc,itime,2,il_top)
vbot=MR_SndVars_metP(iloc,itime,5,il_bot)
vtop=MR_SndVars_metP(iloc,itime,5,il_top)
MR_SndVars_metP(iloc,itime,5,il) = vbot + (vtop-vbot)*(zhere-zbot)/(ztop-zbot)
endif
endif
enddo
MR_SndVars_metP(iloc,itime,5,:) = MR_SndVars_metP(iloc,itime,5,:) + 273.0_sp ! convert C to K
! Wind velocity: do error checking
do il = 1,ulev
! Set NaN's below the surface to zero
if(pres_Snd_tmp(il).gt.SurfPres.or.il.eq.1)then
if(WindVelocity(il).lt.-9998.0_sp) &
WindVelocity(il) = 0.0_sp
else
! For NaN's above the surface, interpolate to the next valid value
if(WindVelocity(il).lt.-9998.0_sp)then
il_bot=il-1
do ill=il,ulev
if(WindVelocity(ill).gt.-9998.0_sp)then
il_top = ill
exit
endif
enddo
zhere=MR_SndVars_metP(iloc,itime,2,il)
zbot=MR_SndVars_metP(iloc,itime,2,il_bot)
ztop=MR_SndVars_metP(iloc,itime,2,il_top)