-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwindgramtt.ncl
2910 lines (2258 loc) · 101 KB
/
windgramtt.ncl
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
; TJ Olney's - Tonino Tarsi 2019 windgram local lapserate and lcl generator
; Requires ncl 5.1.1 for addfiles
; Replace your domains regions sites wlats wlons and wgrids for those below
; can invoked without args and will print the default windgram to x11 display
;
; USAGE: ncl windgramtt.ncl DOMAIN=\"domainname\" [ptop=nn] [rhcut = nn ] [show_[xxx]=[0|1]]
;
; rhcut is a number 1-99 for relative humidity. Default is 95 so >95%rh gets cross hatched.
; Choose your favorite number to cross hatch the time/heights with high enough rh for clouds
; if invoked with ncl DOMAIN=\"domainname\" it will do all the sites in that domain
; ptop is the number of levels, don't go over 3 less than the number of levels in the model, RASP usually 52-3=49
; as is, it calculates lcl three different ways and you must change it in the gsn_csm_addpolymarkers = command if you want to try another
; it has not been completely "error trapped" but it catches some.
;
; New 7/4/2009 posting no blipspots used to gather data, they are recalculated using DrJack functions.
; cldfra_plot added to plot conditionally if there are supposed to be clouds, this is quite experimental, don't know what it will look like in practice.
; hgldj_p used to plot little overturned crescent moons that happen to look a bit like paragliders (couldn't find any better wings)
; using drjack's lcl measure... quite close to lcl1 from lclvl(). You can plot both simply by uncommenting the line
; define default data source for testing purposes
; add command line cldfra_cut decimal for % at which to say it is cloudy enough to try to plot clouds.
;
; Coordinates space fix for windbarbs added. 2011_08_05
; Andrea Barcellona's fix to correct for the map projection and put windbarbs on a
; N-S vertical E-W horizontal coordinate system when the grid is tilted from due N-S
BASEDIR=getenv("BASEDIR")
NCARG_ROOT=getenv("NCARG_ROOT")
if (.not.isvar("DOMAIN")) then
print("Usage: windgramtt.ncl DOMAIN=domainname SITEDATA=sitedata.ncl [ptop=nn] [rhcut = nn ] [show_[xxx]=[0|1]]")
exit
end if
if (.not.isvar("debug")) then ;debug surrounds numerous diagnostic print statements set to anything other than 0 for verbose printing
debug=0
end if
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
external NCL_JACK_FORTRAN "$BASEDIR/GM/LIB/ncl_jack_fortran.so" ;
dayNames = (/"Domenica","Lunedi","Martedi","Mercoledi","Giovedi","Venerdi","Sabato"/)
;--------------------------- function for reading data from a blipmap .dat file. Tonino Tarsi 2019 ---------------
; undef("getblipmap_data")
; function getblipmap_data(bregion:string,bparam:string ,bnumtimes:integer, bnumy:integer,bnumx:integer)
; local filecommand1,fs,a,i,retData
; begin
; retData = new( (/bnumtimes,bnumy,bnumx/), float) ; has a time dimension
; filecommand1= "ls $BASEDIR/"+bregion+"/OUT/" + bparam + ".curr*.data " ;
; fs = systemfunc( filecommand1 ) ;
; do i=0,bnumtimes-1
; print("getblipmap_data - Reading : " + fs(i))
; ;print(dimsizes(retData))
; a = readAsciiTable(fs(i),bnumx,"float",4 )
; ;print(dimsizes(a))
; retData(i,:,:) = a
; end do
; return (retData)
; end
;;; ROUTINE TO SPLIT STRING INTO ARRAY BASED ON INPUT DELIMITER CHARACTER (single character only)
undef("split_string_into_array")
function split_string_into_array( inputstring:string, delimiter:string )
local string_chars, temp_strings, string_array, string_max, nstrings, ichar,ichar1,ichar2
begin
;;; convert delimiter-delimited environmental string into parameter array (what a mess!)
;4testprint: print ( "INPUTSTRING= "+inputstring )
string_chars = stringtochar( inputstring )
string_max = dimsizes(string_chars)-2
;;; temporary array determines maximum parameter names
temp_strings = new( (/300/), string)
ichar1=0
nstrings=-1
do ichar=1,string_max
; allow for either delimiter or no-delimiter at end of string
if( string_chars(ichar).eq.delimiter .or. ichar.eq.string_max) then
nstrings = nstrings +1
if( ichar.eq.string_max) then
ichar2 = string_max
else
ichar2 = ichar-1
end if
temp_strings(nstrings) = chartostring( string_chars(ichar1:ichar2) )
;4testprint: print ( "NAME= "+nstrings+" : "+ichar1+"-"+ichar2+" => "+temp_strings(nstrings) )
ichar1= ichar+1
end if
end do
string_array = new( (nstrings+1), string)
do iname=0,nstrings
string_array(iname) = temp_strings(iname)
end do
return(string_array)
end
;;;---------------------------------------------------------------------------------
operator = "TT"
default_type = "png"
default_ptop = 40 ; usually 20, I changed for TS upper instability
default_rhcut = 95
tilted_grid=0 ; set to 1 if your grid is tilted from NS to cover an oblique area like England or Italy.
one_site = False ; ??? why needed? see near line 369
testndx=0 ; in case a site requested out of range
; Color choice (must be in the palette being used or the closest one will be used)
foreground_color= "black"
background_color= "white"
;background_color= (/.5, .5, .9 /) ; color could be set even with RGB triplet (this example is a deep purple-blue)
label_color= foreground_color
wstar_color= foreground_color
line_color= 10
templine_color= 14
colorpalette = ""
; options to display or not on the windgram
; the check for existing variable allows them to be passed on the command line.
if (.not.isvar("white_clouds")) then
white_clouds=1 ;
end if
if (.not.isvar("use_only_hours")) then
use_only_hours=0 ;
end if
if (.not.isvar("tmp_units")) then
tmp_units = 0 ; change to 0 for C for Celsius, 1 for Fahrenheit; only affects the show_temp variable.
end if
if (.not.isvar("show_cutop")) then
show_cutop = 1 ; CUTOP (top of convection level from cumulus par)
end if
if (.not.isvar("show_sun")) then
show_sun = 1 ; "Solar radiation % from Blipmap
end if
if (.not.isvar("show_hct")) then
show_hct = 0 ; "Cloud top height from TEMF PBL" Requite TEMF schema .. expermental
end if
if (.not.isvar("show_od")) then
show_od = 1
end if
if (.not.isvar("show_blcloudpct")) then
show_blcloudpct = 1
end if
if (.not.isvar("show_cape")) then
show_cape = 1
end if
if (.not.isvar("show_cape_3d")) then
show_cape_3d = 0
end if
if (.not.isvar("show_condense")) then
show_condense = 0
end if
if (.not.isvar("show_cloudfrac")) then
show_cloudfrac = 1
end if
if (.not.isvar("show_wind")) then
show_wind = 1 ; plot the windbarbs
end if
if (.not.isvar("color_wind_barbs")) then
if (show_wind.eq.1) then
color_wind_barbs=0 ; change this if you want the wind barbs to be colored by a scalar of their speed note units
else
color_wind_barbs=1
end if
end if
if (.not.isvar("show_rain")) then
show_rain=1 ; change to 1 if you want to show rain. TJ 6/1/2011
show_this_rain=0
end if
if (.not.isvar("show_lapse")) then
show_lapse = 1 ; local level to level lapse rate
end if
if (.not.isvar("show_wing")) then
show_wing = 1 ; hg/pg 225 ft/min (1.2 m/s) soaring level markers
end if
if (.not.isvar("show_lcl")) then
show_lcl = 1 ; little clouds for lifted condensation level
end if
if (.not.isvar("show_snow")) then
show_snow = 0 ; lowest freezing level
end if
if (.not.isvar("show_pbl")) then
show_pbl = 1 ; boundary layer line
end if
; if ((show_wing.ne.0).or.(show_lcl.ne.0).or.(show_snow.ne.0).or.(show_pbl.ne.0)) ; force lapse if any of the others
; show_lapse = 0
; end if
if (.not.isvar("hide_lapse_colors")) then
hide_lapse_colors=0 ; any other value will make all colors transparent for lapse so it doesn't show but still provides scaffold for other stuff.
end if
if (.not.isvar("show_rh")) then
show_rh = 1 ; relative humidity greater than rh_cut, cheap trick for clouds in the BL; TJ use 95%
end if
if (.not.isvar("show_temp")) then
show_temp = 1 ; temperature contours
end if
if (.not.isvar("show_up_vel")) then
show_up_vel = 1 ; show vertical velocity W* at top of time column set to 0 if no access to DRJack Fortran modules.
end if
if (.not.isvar("show_z_contour")) then
show_z_contour = 0 ; draw a contour line for the altitude attaching pressure level in the morning to corresponding altitude in evening
z_stride=4 ; extra spacing between altitude contour lines
end if
if (.not.isvar("show_p_contour") .and. show_z_contour.ne.0 ) then
show_p_contour= 1 ; draw a contour line for the pressure level in the morning to corresponding altitude in evening
p_stride=z_stride ; extra spacing between altitude contour lines
else
show_p_contour= 0
end if
;; where to find the wrfout files to use as raw data.
;; don't complain about fallbacks when parameters are out of range.
err = NhlGetErrorObjectId()
setvalues err
"errLevel" : "Fatal" ; only report Fatal errors
end setvalues
if (.not.isvar("type")) then
type = default_type
end if
if (.not.isvar("ptop")) then
ptop = default_ptop
end if
if (.not.isvar("rhcut")) then
rhcut = default_rhcut
end if
; this only works if domain names have WINDOW or not
if (isvar("DOMAIN"))
if (isStrSubset(DOMAIN, "WINDOW"))
grid="w2"
else
grid="d2"
end if
end if
; generate house keeping data
offsetstring= systemfunc("date +%z")
pstpdt= systemfunc("date +%Z")
mm_dd=systemfunc("date -u +%b-%d")
fcstday=systemfunc("date +%d' '%B' '%Y")
yyyy_mm_dd=systemfunc("date -u +%Y-%m-%d")
offset = stringtoint(offsetstring)/100
xaxis_msg =""
if (debug.ne.0)
print("")
print("------------------------")
print("Windgrams initialization")
print("------------------------")
print("offsetstring: "+offsetstring)
print("pstpdt: "+pstpdt)
print("mm_dd: "+mm_dd)
print("fcstday: "+fcstday)
print("yyyy_mm_dd: "+yyyy_mm_dd)
print("offset: "+offset)
end if
; Site(s) setup
; domaindata.ncl contains a simple list of domains one per line to be processed and associated with the sites in sitedata.ncl
; PNW-WINDOW
; sitedata.ncl contains long, lat, grid domain and name information each site on a row comma delimted.
; PNW-WINDOW,Iron,w2, -121.94, 48.48,
; must be careful not to let any trailing spaces into sitedata.ncl
; it is read and processed and sorted by region and NS into different arrays by readsitelist.ncl
; I have in the same file a debug mode that generates lists for javascripts and html menus from all the sites to ease web maintenance
if (debug.ne.0)
print("")
print("----------------------------------")
print("Start of readsitelist.ncl printout")
print("----------------------------------")
end if
;load "readsitelist.ncl"
; _____________________read site list__________________________________________________
sitefile="sitedata.ncl"
if (isvar("SITEDATA"))
sitefile = SITEDATA
end if
print("Using " + sitefile)
domains = (/ DOMAIN /)
tmpallsite_params=asciiread(sitefile, -1, "string")
sqsort(tmpallsite_params)
allsite_params=tmpallsite_params
delete(tmpallsite_params)
numdomains=dimsizes(domains) ;number of domains is number of lines in the domain file.
numsites=dimsizes(allsite_params) ; We now have the total number of sites in the database.
print("There are "+numdomains+" domains and a total of "+numsites+" sites in the "+sitefile+" site file.")
site_params=new((/6,numsites /),string) ; ncl array of the original data plus an index number in zero position
sitedomain_ndx=new(numsites,integer) ; should have the index number corresponding to the domain from the domains array
j = 0
do i=0,numsites-1
print(allsite_params(i))
print(str_index_of_substr(allsite_params(i),"#",0))
if ( ismissing(str_index_of_substr(allsite_params(i),"#",0) ) ) then
site_params(0,j)=j ; assign a site number and put it in position zero
site_params(1:,j)=split_string_into_array(allsite_params(i),"," )
site_params(1,j)= DOMAIN
j= j+1
end if
end do
numsites = j
site_params := site_params(:,:j-1)
site_ndx = stringtointeger(site_params(0,:)) ; make it an integer array for use as a subscript
sitedomain = site_params(1,:)
sitename = site_params(2,:)
sitegrid = site_params(3,:)
sitelon = stringtofloat( site_params(4,:) )
sitelat = stringtofloat( site_params(5,:) )
domain=sitedomain
sites=sitename
wlats=sitelat
wlons=sitelon
wgrids=sitegrid
regions=sitedomain
sites_in_domain=new(numdomains,integer)
sites_in_domain(0)=numsites
;______________________________________________________________________________________
if (debug.ne.0)
print("--------------------------------")
print("End of readsitelist.ncl printout")
print("--------------------------------")
end if
;delete(SITES)
;delete(locations) ; not used so delete
numdomains=dimsizes(domains)
numsites=dimsizes(sites)
numregions=dimsizes(regions)
numsitesperdomain=new(numdomains,integer)
wsites=new(numsites,string)
wlon=new(numsites,float)
wlat=new(numsites,float)
if (debug.ne.0)
print("")
print("numdomains: "+numdomains)
print("numsites: "+numsites)
print("numregions: "+numregions)
print("")
end if
do i=0,numdomains-1
numsitesperdomain(i)=num(regions.eq.domains(i)) ; count the sites in each domain
if (debug.ne.0)
print("domain("+i+") = "+domains(i)+" and it has "+numsitesperdomain(i)+" sites held in numsitesperdomain("+i+")" )
end if
end do;
; Discover what we are doing and set up looping parameters for it; what about an "all" case
if (isvar("DOMAIN").and.DOMAIN.eq."all")
if (debug.ne.0)
print ("Doing all Domains ")
end if
numloops=numdomains-1 ; set to loop through domains
else
numloops = 0 ; only using one set of wrfout files
end if
; Begin domains outer loop once for each domain, go through it once if numloops is 0
do nd=0,numloops
if (debug.ne.0)
print("")
print("-----------------------")
print("Start outer domain loop")
print("-----------------------");
print("")
end if
thisdomain_ndx=nd
region=domains(nd) ; regions and domains are the same no window added
if (debug.ne.0)
print("thisdomain_ndx is "+thisdomain_ndx)
print("Top domain: "+region) ; in RASP each named region may have one or two domains
end if ; hi-res grids for domains are designated w2 lo-res d2
; not all but a DOMAIN is being plotted
; code block DOMAIN, extract domain_ndx as pointer
if (isvar("DOMAIN").and.(DOMAIN.ne."all")) ; if domain other than all is passed, find it and set region to it
if (debug.ne.0)
print ("A domain has been passed on the command line: "+DOMAIN+", one of "+numdomains+" possible domains")
end if
do i=0,numdomains-1
if ( (DOMAIN .eq. domains(i)).or.(DOMAIN.eq.str_upper(domains(i))) ) then
region=domains(i)
thisdomain_ndx=i
if (debug.ne.0)
print("The domain "+DOMAIN+" was asked for. Its index is now set to "+thisdomain_ndx+" or thisdomain_ndx")
end if
break
end if
j=i+1
if (j.eq.numdomains)
print("There is no domain known by the name of "+DOMAIN+". Did you mean "+DOMAIN+"-WINDOW perhaps ?")
exit
end if
end do;
end if ; end if (isvar("DOMAIN").and.(DOMAIN.ne."all")) ; if domain other than all is passed, find it and set region to it.
; then do the same for the region
; if a single SITE is passed, find the appropriate members of the arrays
; code block SITE
if (isvar("SITE"))
do i=0,numsites-1 ; go thru sites until you find the one passed.
if ( (SITE .eq. sites(i)).or.(SITE.eq.str_upper(sites(i))) ) then
wsites(0)=sites(i)
wlon(0)=wlons(i)
wlat(0)=wlats(i)
wgrid=wgrids(i)
region=regions(i)
DOMAIN=regions(i)
thissite_ndx=site_ndx(i)
one_site= True
numtodo=1
loopstart=thissite_ndx
print("Code block SITE: single site passed "+SITE+" which is site "+sites(site_ndx)+", value of thissite_ndx is "+thissite_ndx)
break ; found it nothing to gain by carrying on.
end if
end do;
testndx=i+1
else
one_site= False
end if ; end isvarsite
if (testndx.ge.numsites)
print("You have entered a request for a site, "+SITE+" not in the sitedata.ncl file. ")
print(sites)
exit
end if
; nothing passed, so just do the default to the workstation
; code block NOTHING PASSED DEFAULT WINDGRAM
if (.not.isvar("DOMAIN").and.(.not.isvar("SITE") )) then
print ("no SITE or DOMAIN variable using defaults to x11 workstation. Default parameters")
thissite_ndx=34 ; choose the default index set all else from that.
wsites(0)=sites(thissite_ndx)
wlon(0)=wlons(thissite_ndx)
wlat(0)=wlats(thissite_ndx)
wgrid=wgrids(thissite_ndx)
region=regions(thissite_ndx)
DOMAIN=regions(thissite_ndx)
one_site= True
numtodo=1
loopstart=thissite_ndx
SITE=sites(thissite_ndx)
if (type.eq."ncgm")
type="x11"
end if
print("")
print("---------------- Running default windgram ----------------");
print("DOMAIN = "+DOMAIN)
print("SITE = "+SITE+" whose index is "+thissite_ndx)
print("")
; something was passed to the script so use it to decide what to plot
; if DOMAIN, region has already been set if SITE parmams have been set
; so remaining case is DOMAIN without SITE so arrays of sites in the domain must be set
; DOMAIN is region(domain_ndx) is also region
else
print("There is a domain = "+DOMAIN )
; no domain so see if there is a site
; they are the whole possible arrays but re assign them by checking whether site is in domain
if (one_site.eq.False) then
testdomain= DOMAIN
; Get new arrays for just this domain/region - Why ???
jay=0
do i=0, numsites-1 ; must check them all to find which belong
if (regions(i).eq.region) then ; .or.(regions(i)+"-WINDOW".eq.region)) - WINDOW is now part of the domain/region name.
print("for sitename(i) "+sitename(i)+" site_ndx "+site_ndx(i)+" sitedomain(i) "+sitedomain(i)+" matches "+testdomain+" ")
wsites(jay)=sites(i)
print ("sites("+i+")="+sites(i) )
wlon(jay)=wlons(i)
wlat(jay)=wlats(i)
wgrid=wgrids(i)
if (debug.ne.0)
print ("sites("+i+")="+sites(i) )
print( "wlon(jay)="+wlon(jay)+" wlons(i) = "+wlons(i)+"wlat(jay)="+wlat(jay)+" wlats(i) = "+wlats(i) )
end if
if (jay.eq.0)
thissite_ndx=i
end if
jay=jay+1
end if
if (jay.eq.1)
one_site=True
end if
end do;
end if
end if
; end if (.not.isvar("DOMAIN").and.(one_site.eq.False ))
if (.not.isvar("numtodo"))
numtodo=sites_in_domain(thisdomain_ndx)
end if
if (numtodo.eq.1)
one_site=True
end if
if (debug.ne.0)
print("")
;print("Generating windgrams for region "+region+", grid "+wgrid+" has "+numtodo+" sites")
;print("domain_ndx = "+domain_ndx)
print("")
print("wsites,wlon,wlat,wgrid,region")
do i=0,numtodo-1
print(wsites(i)+","+wlon(i)+","+wlat(i)+","+wgrid+","+region)
end do;
end if
if (debug.ne.0)
print(" ")
print("Plotting parameters: ")
print(" ")
print("tilted_grid = "+tilted_grid)
print("ptop = "+ptop)
print("Output plot type = "+type)
print("show_lapse = "+show_lapse)
print("show_wind = "+show_wind)
print("color_wind_barbs = "+color_wind_barbs)
print("show_rh = "+show_rh)
print("rhcut = "+rhcut)
print("show_temp = "+show_temp)
print("debug = "+debug)
print("tmp_units = "+tmp_units)
print("show_condense = "+show_condense)
print("show_rain = "+show_rain)
print("show_wing = "+show_wing)
print("show_lcl = "+show_lcl)
print("show_snow = "+show_snow)
print("show_pbl = "+show_pbl)
print("show_cloudfrac = "+show_cloudfrac)
print("hide_lapse_colors = "+hide_lapse_colors)
print("show_up_vel = "+show_up_vel)
print("show_z_contour = "+show_z_contour)
print("show_p_contour = "+show_p_contour)
end if
; Data ingest from wrfout files and digest
if ( use_only_hours.eq.1 ) then
filecommand1= "ls $BASEDIR/"+region+"/wrfout_d02_* | egrep -v ':30|:03:|04:|05:|06:|17:|18:'" ; default CUCCO set at top
else
filecommand1= "ls $BASEDIR/"+region+"/wrfout_d02_* | egrep -v ':03:|04:|05:|06:|17:|18:'" ; default CUCCO set at top
end if
print("")
print("WRF data files are found as "+BASEDIR+"/"+region+"/wrfout_d02_*")
print("Forecast date is UTC "+yyyy_mm_dd+" "+offset+" to get "+pstpdt )
print("")
; if ( wgrids(idx) .eq. "w2" ) then ; obsolete
; region=regions(i)+"-WINDOW"
; end if
; set at top
fs = systemfunc( filecommand1 ) ; fs becomes an array of full filenames, no more needs deleting between domains
onefile=fs(0) ; get a single file from which to get geographic data, time labels, etc
print("Processing these data files: ")
print(fs+" ")
print("")
print("Fixed parameters, ter, z, etc, being looked up in the file :")
print(onefile+" ")
a = addfiles(fs+".nc","r") ; addfiles grabs all the files in the list above at once.
timesinfiles = a[:]->Times ; get date & time from wrfout files
b = addfile(onefile+".nc","r") ; addfiles grabs the file used for fixed parameters
p=a[:]->P ; perturbation pressure (Pa) for a time,ptop,long,lat array
pb=a[:]->PB ; base state pressure (Pa) for a time,ptop,long,lat array
press=p ; bring along metadata
press = p + pb ; total pressure
pmb=press/100 ; convert total pressure in hPa
pdims=dimsizes(press)
numlevels=pdims(1)
numtimes=pdims(0)
numy=pdims(2)
numx=pdims(3)
taus=ispan(0,numtimes-1,1) ; coordinate variable for Time dimension
ilevels=ispan( 0,numlevels-1,1)
tkall=a[:]->T ; perturbation potential temperature
tkall = tkall + 300 ; as per wrfuser
tk = wrf_tk( press , tkall ) ; temperature Kelvin
tc=tk
tc = tk - 273.16
tc@units = "C"
ter=a[:]->HGT ; meters above sea level
pblh=a[:]->PBLH ; meters pbl thickness
pblt= ter + pblh ; height of top of pbl meters
tmpu =a[:]->U
tmpv =a[:]->V
u = wrf_user_unstagger(tmpu,tmpu@stagger)
v = wrf_user_unstagger(tmpv,tmpv@stagger)
QVAPOR = a[:]->QVAPOR
td=wrf_td(press,QVAPOR) ; get the dewpoint everywhere, use for all sites
rh = wrf_rh(QVAPOR, press, tk) ; get the relative humidity everywhere, use for all sites
; PH = a[:]->PH
; PHB = a[:]->PHB
; PH = PH + PHB
; PSFC = a[:]->PSFC ; don't use for RASP time 0 is always = 0
; z = wrf_user_unstagger(PH,PH@stagger)/9.81 ; On full (w) levels
z = wrf_user_getvar(a,"z",-1) ; ToninoTarsi 2019 Get z direcly from wrf in meters . http://www.ncl.ucar.edu/Document/Functions/WRF_arw/wrf_user_getvar.shtml
dimv = dimsizes(z)
zmeter = z ; used below in meters
z@units = "m"
z!0 = "Time"
z!1 = "levels"
z!2 = "south_north"
z!3 = "west_east"
z&Time = taus
zmeter@units = "m"
zmeter!0 = "Time"
zmeter!1 = "levels"
zmeter!2 = "south_north"
zmeter!3 = "west_east"
zmeter&Time = taus
mslp = wrf_slp( zmeter, tk, press, QVAPOR )
;mslp = wrf_user_getvar(a,"slp",-1) ; ToninoTarsi 2019
if ( show_hct.ne.0) then ; ToninoTarsi 2019 check if we run TEMF schema
if ( isfilevar(b,"HCT_TEMF") ) then
hct = a[:]->HCT_TEMF
else
show_hct = 0
end if
end if
if ( show_cutop.ne.0) then ; ToninoTarsi 2019 check if CUTOP present in WRF data
if ( isfilevar(b,"CUTOP") ) then
cutop = a[:]->CUTOP
cubot = a[:]->CUBOT
else
print("Warning ... CUTOP not founf in WRF output")
show_cutop = 0
end if
end if
dims = dimsizes(timesinfiles)
times = new(dims(0),string) ; string representation of the times in the files
caption = new(dims(0),integer)
lst = new(dims(0),integer) ; put local time into this for labels
lst_string = new(dims(0),string) ; put local time in string variable for labels
lst_string_striped = new(dims(0),string) ; put local time in string variable for labels
lst_is_hour = new(dims(0),integer)
start_time_margin = 1.0
if ( use_only_hours.eq.1 ) then
start_time_margin = 0.5
end if
do i=0,dims(0)-1
times(i) = chartostring(timesinfiles(i,8:15))
tmp = chartostring(timesinfiles(i,11:12))
caption(i) = stringtoint(tmp) ; +"."+chartostring(t(2,10:12))
lst(i) = caption(i)+offset ; convert UTC to local standard time
; make times 12 hour local clock
if ( lst(i).lt.1 ) ; negative numbers add 24
lst(i)= lst(i)+24
end if
; add ".00" to local time
lst_string(i) = sprinti("%0.2i", lst(i)) + ":" + chartostring(timesinfiles(i,14:15))
lst_string_striped(i) = sprinti("%0.2i", lst(i)) + chartostring(timesinfiles(i,14:15))
lst_is_hour(i) = str_index_of_substr(lst_string(i),":00",-1)
ora_solare = True
if ( ora_solare ) then
if ( lst_string(i).eq."09:00") then
start_time_index = i - start_time_margin
end if
if ( lst_string(i).eq."17:00") then
stop_time_index = i + start_time_margin
end if
else
if ( lst_string(i).eq."10:00") then
start_time_index = i - start_time_margin
end if
if ( lst_string(i).eq."18:00") then
stop_time_index = i + start_time_margin
end if
end if
end do
undef("getblipmap_data")
function getblipmap_data(bregion:string,bparam:string ,bnumtimes:integer, bnumy:integer,bnumx:integer)
local filecommand1,fs,a,i,retData
begin
retData = new( (/bnumtimes,bnumy,bnumx/), float) ; has a time dimension
do i=0,bnumtimes-1
filecommand1 = BASEDIR+bregion+"/OUT/" + bparam + ".curr." + lst_string_striped(i) + "lst.d2.data"
;print(filecommand1)
if ( fileexists(filecommand1) ) then
a = readAsciiTable(filecommand1,bnumx,"float",4 )
retData(i,:,:) = a
else
retData(i,:,:) = default_fillvalue("float")
;print("no data for " + bparam + " at " + lst_string_striped(i))
end if
end do
return (retData)
end
if (show_rain.ne.0) then
rain1_multitime = getblipmap_data(region,"rain1",numtimes, numy,numx)
end if
if (show_sun.ne.0) then
sfcsunpct_multitime = getblipmap_data(region,"sfcsunpct",numtimes, numy,numx)
end if
if (show_cape_3d.ne.0) then
c3d = wrf_user_getvar(a,"cape_3d",-1) ; ToninoTarsi 2019
cape3d_multitime = c3d(0,:,:,:,:)
cin3d_multitime = c3d(1,:,:,:,:)
cape2d_multitime = dim_max_n(cape3d_multitime,1)
cin2d_multitime = dim_max_n(cin3d_multitime,1)
cinfo = wrf_user_getvar(a,"cape_2d",-1)
mcape = cinfo(0,:,:,:)
mcin = cinfo(1,:,:,:)
lcl = cinfo(2,:,:,:)
lfc = cinfo(3,:,:,:)
end if
if (show_cape.ne.0) then
cape_multitime = getblipmap_data(region,"cape",numtimes, numy,numx)
end if
if (show_blcloudpct.ne.0) then
blcloudpct_multitime = getblipmap_data(region,"blcloudpct",numtimes, numy,numx)
end if
if ( show_od.ne.0) then
zblcldif_multitime = getblipmap_data(region,"zblcldif",numtimes, numy,numx)
zblcl_multitime = getblipmap_data(region,"zblcl",numtimes, numy,numx)
end if
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; To calculate Wstar up velocity
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if ( .not. isvar("vhf") ) then
if ( .not. isvar("hfx") ) then
hfx = a[:]->HFX ;,time) ; for sfc. sensible heat flux in w/m2
end if
; convert latent heat flux into additional virtual temperature heat flux
; 0.61*cp/L=0.61*(1006J/Kkg)/(2.502x106J/kg)=0.000245268
vhf = a[:]->LH
; DRJACK routines calculate for only one time period so vhf is only 2d, not 3d
; must calculate for each period
; hcrit wstar sfclcl directly
do i=0,numtimes-1
NCL_JACK_FORTRAN::minlimit2d( vhf(i,:,:), 0.0, numx,numy )
vhf(i,:,:) = hfx( i,:,:) + 0.000245268*(tc(i,0,:,:)+273.16)*vhf(i,:,:)
end do;
end if ; That should fully populate vhf
; s_n w_e
wstar = new( (/numtimes,numy,numx/), float) ; has a time dimension
wstar_1 =new( (/numy,numx/), float) ; no time dimension, returned by NCL_JACK_FORTRAN
hcrit = new( (/numtimes,numy,numx/), float) ; has a time dimension
hcrit_1 =new( (/numy,numx/), float) ; no time dimension, returned by NCL_JACK_FORTRAN
sfclclheight = new( (/numtimes,numy, numx/), float)
sfclclheight_1 = new( (/numy, numx/), float)
hglider_1 = new( (/numy, numx/), float)
hglider = new( (/numtimes,numy, numx/), float)
; NCL_JACK_FORTRAN :: calc_wstar( vhf,pblh, isize,jsize,ksize, wstar )
; use DrJack's 2d functions for other parameters of interest -- loop over numtimes
do i=0,numtimes-1
if (show_up_vel.ne.0.or.show_wing.ne.0) then
NCL_JACK_FORTRAN :: calc_wstar( vhf(i,:,:),pblh(i,:,:), numx,numy,numlevels, wstar_1 )
wstar(i,:,:)=wstar_1 ; put result for this iteration into the i time slot for all wstar
end if
if (show_wing.ne.0) then
NCL_JACK_FORTRAN :: calc_hcrit( wstar_1, ter(i,:,:), pblh(i,:,:), numx,numy, hcrit_1 )
hcrit(i,:,:)=hcrit_1
; hcrit(i,:,:)=hcrit_1*3.28084 ; to feet ; put result for this iteration into the i time slot
end if
if (show_lcl.ne.0 .or. show_wing.ne.0) then ; need it to calculate hcrit of wing even if lcl not shown
NCL_JACK_FORTRAN :: calc_sfclclheight( press(i,:,:,:), tc(i,:,:,:), td(i,:,:,:), zmeter(i,:,:,:), ter(i,:,:), pblh(i,:,:), numx,numy,numlevels, sfclclheight_1 )
sfclclheight(i,:,:)= sfclclheight_1 ; put result for this iteration into the i time slot
end if
end do;
; sfclclheight=sfclclheight*3.28084 ; meters to feet
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Wstar calculated for whole region
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
llres = True
llres@ReturnInt = True ; return integer values
llres@MAP_PROJ =b@MAP_PROJ
llres@TRUELAT1 =b@TRUELAT1
llres@TRUELAT2 =b@TRUELAT2
llres@STAND_LON =b@STAND_LON
llres@DX = b@DX
llres@DY = b@DY
XLAT = b->XLAT(0,:,:)
XLONG = b->XLONG(0,:,:)
day=chartostring(timesinfiles(1,0:9))
plotday=systemfunc("date --date="+day+" +%d' '%b' '%Y")
print ("Plotday is "+plotday)
print ("For day "+day+" for "+region+" numtimes is number of times available = "+numtimes) ; leave in report
theYear = stringtointeger(chartostring(timesinfiles(1,0:3)))
theMonth = stringtointeger(chartostring(timesinfiles(1,5:6)))
theDay = stringtointeger(chartostring(timesinfiles(1,8:9)))
theWeekNumber = day_of_week(theYear,theMonth,theDay)
theDayName = dayNames(theWeekNumber)
feet=new(pdims(1),float)
; End of whole domain definitions
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Setup sites loop pull out bits for the site location
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; debug start ;;;
if (debug.ne.0)
print("top of loop loopstart = "+thissite_ndx+" index of first site in domain after sorting sites by domain, with "+numtodo+" sites in the region to do \n\n" )
end if
;;; debug end ;;;
if (numtodo.eq.1) then
loopstart=thissite_ndx
;;; debug start ;;;
if (debug.ne.0)
print(" before being used as subscript, loopstart = "+loopstart)
end if
;;; debug end ;;;
loopend=loopstart ; only for onesite case
wlon(loopstart)=wlons(loopstart)
wlat(loopstart)=wlats(loopstart)
wsites(loopstart)=sites(loopstart)
print ("Running a one_site case for site "+thissite_ndx+", "+wsites(loopstart) )
else
loopstart=thissite_ndx
loopend=loopstart+numtodo-1
print ("Running a domain case for "+region)
end if
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Actual start of do loop for sites in one region
; ns tracks which site
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do ns=loopstart,loopend ; for one loxY locX at a time.
outdir=BASEDIR+"/"+region+"/OUT/" ; windgram output directory,
wpt_code := str_split(sites(ns), " ")
if ( dimsizes(wpt_code).eq.0 ) then
continue
end if
windgram=outdir+wpt_code(0)
;windgram=outdir+day+"_"+sites(ns)
;;; debug start ;;;
if (debug.ne.0) then
print(" ")
print(" ")
print(" loopstart = "+loopstart)
print("--------------------top of loop---------------------")
print(" loopend = "+loopend)
print(" ")
print(" ")
print("windgram fullpathname = "+windgram+" output to "+type)
end if
;;; debug end ;;;
; Set up the Workstation for this site
; must be destroyed at end of loop must be within the loop color array can be outside though
; wks = gsn_open_wks(type , windgram) ; open a workstation as type with constructed name
type@wkWidth = 800