File: mv_vapor_prep.mv

package info (click to toggle)
metview 5.10.2-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 242,296 kB
  • sloc: cpp: 437,117; ansic: 41,433; xml: 19,944; f90: 13,059; sh: 6,562; python: 3,953; yacc: 1,774; lex: 1,121; perl: 701; makefile: 92
file content (1024 lines) | stat: -rw-r--r-- 26,833 bytes parent folder | download | duplicates (5)
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
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
#!/usr/bin/magml
# Metview Macro

#  **************************** LICENSE START ***********************************
# 
#  Copyright 2018 ECMWF. This software is distributed under the terms
#  of the Apache License version 2.0. In applying this license, ECMWF does not
#  waive the privileges and immunities granted to it by virtue of its status as
#  an Intergovernmental Organization or submit itself to any jurisdiction.
# 
#  ***************************** LICENSE END ************************************
# 
#---------------------------------------
# Define global parameters
#---------------------------------------

global VDF_FILE=""
global OUTPUT_DIR=""
global WORK_DIR=""
global LOG_FILE=""

global STEP_NUMBER=-1
global SURF_PARAMS=""
global UPPER_PARAMS=""
global VERTICAL_GRID=""
global ELEVATION_PARAM=""
global MIN_VERTICAL="0"
global MAX_VERTICAL="0"
global VDF_LEVEL=""

global AREA=""
global GRID=""

global INPUT_GRB=nil
global DATE_LST=nil
global LEV_LST=nil
global LEV_TYPE=""
global LON_WEST=0
global LON_EAST=0
global LAT_NORTH=0
global LAT_SOUTH=0
global NX=0
global NY=0
global NZ=0
global REVERT_LAT=0
global REVERT_LON=0
global errH=0
global infoH=0

#--------------------------------
# The vapor setup command
#--------------------------------

global VAPOR_SETUP_CMD=". $MV_VAPOR_BIN/vapor-setup.sh"

#----------------------------------------------------
#The exe to convert grib into raw binary for vapor
#----------------------------------------------------

global TORAW_EXE=getenv("METVIEW_BIN") & "/gribToVapor"

#---------------------------------------------------------------
# Read arguments
#---------------------------------------------------------------

#Get arguments
global ARGS
ARGS=arguments()
global CURRENT_ARG_INDEX=1

if count(ARGS) = 0 then
	log_fail("No arguments specified!")
end if

print("Arguments: ",ARGS)

#Get output dir 
OUTPUT_DIR=ARGS[1]
WORK_DIR=OUTPUT_DIR

#Get vdf filename
VDF_FILE=nextArg() & ".vdf"

#Get vdf check mode
val=nextArg()
doCheck=0
if val = "1" then
    doCheck=1
    if exist(VDF_FILE) then
    	log_info("VDF file already exists! No conversion is performed!")
    	stop()
    end if	
end if 

#Get vdf refinement level
VDF_LEVEL=nextArg()

#Get vertical grid params
VERTICAL_GRID=nextArg()
if VERTICAL_GRID = "LAYERED" then
     ELEVATION_PARAM=nextArg()
	 MIN_VERTICAL=nextArg()
	 MAX_VERTICAL=nextArg()
end if

#Get interpoltaion params
val=nextArg()
doInterpolate=0
if val = "1" then
    doInterpolate=1
    AREA=parse(nextArg(),"/")
    GRID=parse(nextArg(),"/")
end if 

#Get step number
STEP_NUMBER=number(nextArg())

#Get parameter list
SURF_PARAMS=parse(nextArg(),"/")
#Surface params are optional
if SURF_PARAMS[1] = "_UNDEF_" then
	 SURF_PARAMS = nil
end if

UPPER_PARAMS=parse(nextArg(),"/") 

#Get working dir
WORK_DIR=nextArg()

#Get input files
inputMode=nextArg()
if inputMode = "ICON" then
    INPUT_GRB=nil
	while hasNextArg() do 
			INPUT_GRB=INPUT_GRB & read(nextArg())
	end while		
else
	log_fail("InputMode=" & inputMode & " not implemented!")
end if	

#---------------------------------------
#Log input parameters
#---------------------------------------

log_info("Input arguments --->")
log_info(tab & "VDF_FILE=" & VDF_FILE)
log_info(tab & "OUTPUT_DIR=" & OUTPUT_DIR)
log_info(tab & "SURF_PARAMS=" &  listToString(SURF_PARAMS))
log_info(tab & "UPPER_PARAMS=" & listToString(UPPER_PARAMS))
log_info(tab & "VERTICAL_GRID=" & VERTICAL_GRID)
if VERTICAL_GRID = "LAYERED" then
     log_info(tab & "ELEVATION_PARAM=" & ELEVATION_PARAM)
	 log_info(tab & "MIN_VERTICAL=" & MIN_VERTICAL)
	 log_info(tab & "MAX_VERTICAL=" & MAX_VERTICAL)
end if
if doInterpolate = 1 then
    log_info(tab & "Interpolate to AREA=" & listToString(AREA) & " GRID=" & listToString(GRID))
end if
log_info(tab & "STEP_NUMBER=" & STEP_NUMBER)
log_info(tab & "VDF_LEVEL=" & VDF_LEVEL)

#---------------------------------------
# Check directories
#---------------------------------------

if not(exist(OUTPUT_DIR)) then
   if shell("mkdir -p " & OUTPUT_DIR) <> 0 then
      print("Could not create output directory!")
      print("Shell command failed: mkdir -p " & OUTPUT_DIR)
      fail()
   end if   
end if

if not(exist(WORK_DIR)) then
   if shell("mkdir -p " & WORK_DIR) <> 0 then
      print("Could not create working directory!")
      print("Shell command failed: mkdir -p " & WORK_DIR)
      fail()
   end if   
end if

#------------------------------------
# Interpolation
#------------------------------------

#Should be done only for the required fields
if doInterpolate = 1 then	 
  	INPUT_GRB=read(
         	data : INPUT_GRB, 
		 	area: AREA,
		 	grid: GRID       
       	)
end if
  
#------------------------------------
# Get vertical structure 
#------------------------------------
getVertical()

#------------------------------------
# Get dates and times
#------------------------------------
getDateTime()

#------------------------------------
# Get geography
#------------------------------------
getGeography()

#------------------------------------
# Sorting
#------------------------------------

if (LEV_TYPE = "ml" and ELEVATION_PARAM <> "119") or LEV_TYPE = "pl" then
   INPUT_GRB=sort(INPUT_GRB,["param","date","time","step","levelist"],["<","<","<","<",">"])
#   #Ig=sort(INPUT_GRB,"parameter")
#   #INPUT_GRB=sort(INPUT_GRB)
else
   INPUT_GRB=sort(INPUT_GRB,["param","date","time","step","levelist"],["<","<","<","<","<"])
end if

#------------------------------------
# Check vertical parameters
#------------------------------------

if VERTICAL_GRID = "LAYERED" then
	prepareElevation()
end if	

#-------------------------
# Create the vdf file
#-------------------------
createVdf()  

#-------------------------
# Process 2D params
#-------------------------

log_info("Processing 2D params --->")

loop v in SURF_PARAMS 
  
  log_info(tab & "param: " & v)
  
  g=read(
         	data : INPUT_GRB, 
		 	levtype: "sfc",
		 	param: v		       
       	)
  if g =nil then
  	   log_fail("No data found for param: " & v) 	   
  end if	
  
  varName=checkVarName(v,"2D") 
  processVar(g,varName,"2D")
      
end loop   

#-------------------------
#Process 3D params
#-------------------------
log_info("Processing 3D params --->")

loop v in UPPER_PARAMS 
  
  log_info(tab & "param: " & v)
  
  g=read(
         data : INPUT_GRB, 
		 levtype: LEV_TYPE,
		 param: v		       
       )
  if g =nil then
  	   log_fail("No data found for param: " & v) 	   
  end if	
        
  varName=checkVarName(v,"3D")      
  processVar(g,varName,"3D")
    
end loop   

stop(0)

#======== Functions ===================================

#=======================================================
# Read the all the validity dates for the first available
# variable. The result is stored in DATE_LST
#=======================================================

function getDateTime()

 log_info("Check dates --->")
 
 if SURF_PARAMS <> nil then
     par=SURF_PARAMS[1]
 else
     par=UPPER_PARAMS[1]
 end if
 
 g=read(data : INPUT_GRB, param: par)
 
 DATE_LST=nil
 loop fg in g	
        
     if STEP_NUMBER = -1 or  count(DATE_LST) <  STEP_NUMBER then
    
        vd=[0,0,0]
     	vdLst=grib_get(fg,["validityDate:l","validityTime:l","step:l"])     	
        if count(vdLst) = 1 then
            vd=vdLst[1]
        end if    
        vdAsDate=valid_date(fg)
     	d=(date: vd[1],time: vd[2], step: vd[3],dd: vdAsDate)        
        
     	if count(DATE_LST) =0 then
     		DATE_LST = DATE_LST & [d]
 	 	else 
 	    	found = 0
 	 		loop item in DATE_LST
 	 	   		if item["date"] = d["date"] and 
 	 	      		item["time"] = d["time"] then	 	     
 	 	          		found = 1 
 	 	   		end if 
 	 		end loop 
 	 	
 	 		if found = 0 then
 	 			DATE_LST = DATE_LST & [d]
 	 		end if        
 	 	end if
 	 end if
 	 		
 end loop		
 
 log_info(tab &"Validity dates:")
 for i=1 to count(DATE_LST) do
 	log_info(tab & i & ") date=" & DATE_LST[i].date & ", time=" & DATE_LST[i].time)
 end for	
 
 #print("dates")
 #print(DATE_LST)
 
end getDateTime

#=======================================================
# Read the all the level related parameters from the first
# upper level parameter.
#=======================================================

function getVertical()

  log_info("Check vertical --->")
  
  par=UPPER_PARAMS[1]
  if par = "lnsp" then
     if count(UPPER_PARAMS) < 2 then
  		log_fail("No upper level params specified!")
  	 end if	
     par=UPPER_PARAMS[2]
  end if
     
  g=read(data : INPUT_GRB, param: par)
  levType="sfc"
  cnt=1
  while levType = "sfc" and cnt <> count(g) do
  	  levType=grib_get_string(g[cnt],"levelType")
  	  cnt=cnt+1
  end while
  
  if levType = "sfc" then
  	 log_fail("No vertical level is found for param: " & par)
  end if	
  
  #Set the level type
  LEV_TYPE=levType
 
  #Generate the level list 
  loop fg in g	
     lev=grib_get_string(fg,"level")
     if count(LEV_LST) = 0 then
     	LEV_LST=LEV_LST & [lev]
     else if lev not in LEV_LST then
        LEV_LST=LEV_LST & [lev] 	
     end if
  end loop
  
  NZ=count(LEV_LST)
   
  log_info(tab & "Level type: " & LEV_TYPE)
  log_info(tab & "Level num: " & NZ)
  levStr=""
  for i=1 to count(LEV_LST) do
  	if i = 1 then
  		levStr=levStr & LEV_LST[i]
  	else
  		levStr=levStr & "," & LEV_LST[i]
  	end if
  end for
  log_info(tab & "Levels: " & levStr)
  
  			
  
end getVertical

function getGeography()

	log_info("Check geography ---> ")
	 
	g=INPUT_GRB[1]
	
	lat1=grib_get_double(g,"latitudeOfFirstGridPointInDegrees")
	lat2=grib_get_double(g,"latitudeOfLastGridPointInDegrees")
	
	if lat1 > lat2 then
		LAT_NORTH=lat1
		LAT_SOUTH=lat2
	else
		LAT_NORTH=lat2
		LAT_SOUTH=lat1	
	end if
		
	lon1=grib_get_double(g,"longitudeOfFirstGridPointInDegrees")
	lon2=grib_get_double(g,"longitudeOfLastGridPointInDegrees")
	
	#log_info(tab & "lon1/lon2: " & lon1 & "/" &lon2)
		
	if lon1 > 180 then
		lon1 = lon1-360
	end if	
	if lon2 > 180 then
		lon2 = lon2-360
	end if
	
	#log_info(tab & "lon1/lon2: " & lon1 & "/" &lon2)

	if lon1 < lon2 then
		LON_WEST=lon1
		LON_EAST=lon2
	else
		LON_WEST=lon2
		LON_EAST=lon1
	end if
		
	NX=grib_get_long(g,"Ni")
	NY=grib_get_long(g,"Nj")
	
	log_info(tab & "nx/ny: " & NX & "/" & NY)
	log_info(tab & "N/W/S/E: " & LAT_NORTH & "/" & LON_WEST & "/" & LAT_SOUTH & "/" & LON_EAST)	
    #log_info(tab & "revert lat: " & REVERT_LAT)
    #log_info(tab & "revert lon: " & REVERT_LON)

end getGeography
	
#====================================================
# Check and modify an original parameter names to 
# fit it to VAPOR's needs
#====================================================		

function checkVarName(varName,varType:string)

  if varName = "z" then  
    if varType = "3D" then
  		return "ELEVATION"
	else if varType = "2D" then
		return "HGT"
	end if
  end if
  	
  if ELEVATION_PARAM = varName then
  	  return "ELEVATION"	
  end if	
  	
  numberLst=["0","1","2","3","4","5","6","7","8","9"]
  if type(varName) = "number" then
  		return "v" & string(varName)
  else		
  	s = substring(varName,1,1)
  	if s in numberLst then
     	return "v" & varName
    end if
  end if  
  
  return varName
end checkVarName     

#====================================================
# Create the VDF file
#====================================================

function createVdf() 
	
  #Write all the timestamps into an ascii file.
  #vdfcreate will need it!	  
  fName=WORK_DIR & "/tsfile.txt"
  f=file(fName) 
  for i=1 to count(DATE_LST) do  
    vdate=DATE_LST[i]["date"]
    vtime=DATE_LST[i]["time"]	   
    write(f,vdate &  vtime, newline)
  end for    
  f=0  
  
  #Generate 2d parameter string for vdfcreate
  par2D=""
  for i=1 to count(SURF_PARAMS) do
	if i=1 then
		par2D=checkVarName(SURF_PARAMS[i],"2D")
	else
		par2D=par2D & ":" & checkVarName(SURF_PARAMS[i],"2D")	
	end if
  end for
	
  #Generate 2d parameter string for vdfcreate
  par3D=""
  for i=1 to count(UPPER_PARAMS) do
	if i=1 then
		par3D=checkVarName(UPPER_PARAMS[i],"3D")
	else
		par3D=par3D & ":" & checkVarName(UPPER_PARAMS[i],"3D")	
	end if
  end for
	
  fLog=WORK_DIR & "/log.txt"
  
  if VERTICAL_GRID = "LAYERED" then
     gridType="layered"
  else
  	 gridType="regular"
  	 if LEV_TYPE = "ml" or LEV_TYPE = "pl" then
  	     MIN_VERTICAL="-" & LEV_LST[1]
  	     MAX_VERTICAL="-" & LEV_LST[count(LEV_LST)]
  	 else
  	     MIN_VERTICAL=LEV_LST[1]
  	     MAX_VERTICAL=LEV_LST[count(LEV_LST)]
  	 end if    
  end if 	    
  
  log_info("---> Create VDF file with vdfcreate!")
  
  #Create VDF file
  geoFactor=111177
  cmd=VAPOR_SETUP_CMD & ";" &
      " cd " & OUTPUT_DIR & ";" & 
      " vdfcreate  -dimension " & NX & "x" & NY & "x" & NZ &
          " -gridtype " & gridType &
          " -mapprojection \"+proj=latlong  +ellps=sphere\" " &
          " -level " & VDF_LEVEL &
          " -numts " & count(DATE_LST) &
 #         " -usertimes tsfile.txt " &
          " -extents " & LON_WEST*geoFactor & ":" &
                LAT_SOUTH*geoFactor & ":" &
               	MIN_VERTICAL & ":" &
               	LON_EAST*geoFactor & ":" &
                LAT_NORTH*geoFactor & ":" &
                MAX_VERTICAL & 
          " -vars3d " & par3D
   
   if par2D <> "" then                
          cmd=cmd & " -vars2dxy " & par2D
   end if	           
          
   cmd=cmd & " " & VDF_FILE &  " > " & fLog            

  print(cmd)
  if shell(cmd) <> 0 then
  	  log_error("vdfcreate failed!")
  	  log_error("command: " & cmd)
  	  #log_error("output " & read(fLog))
  	  fail()
  else
  	log_info("command: " & cmd)		
  end if
  	
    #------------------------------------------------------
    # write timestamps into the vdf file. vdfcreate
    # cannot do this so we write the timestamps directly 
    # into the vdf file.
    #-------------------------------------------------------
 
    #Create a list with all the timestamps
    tsLst=nil
    for i=1 to count(DATE_LST) do          
        dd=DATE_LST[i]["dd"]
        tsLst=tsLst & [string(dd,"yyyy-mm-dd") & " " & string(dd,"HH:MM") &
             " (+" & DATE_LST[i].step & ")"]
    end for
  
    #add the timestamp into the cdf file. The vdf file must have a 
    #  <UserTimeStampString ...>
    #  </UserTimeStampString> 
    #block
  
    log_info("Write timestamps into vdf file")
  
    vdfPath=OUTPUT_DIR & "/" & VDF_FILE
    
    #we copy the vdf file into a tmp file otherwise metview 
    #cannot read it as a simple text file!
    vdfPathTmp=tmpfile() 	  
    shell("cp " & vdfPath & " " & vdfPathTmp)
   
    txt=read(vdfPathTmp)
    if txt = nil then
        log_error("Could not read vdf file: " & vdfPath)
        log_error("Timestamps will not be added to vdf file")
        return
    end if 
  
    modTxt=""
    foundUTSS=0
    tagStart=0
    ts=1
    
    #Iterate through the vdf file line by line. For each time step
    #the vdf file contains a separate block. We try to update the 
    #<UserTimeStampString> tag within each block with our time
    #stamp string.
    for i=1 to count(txt) do
    
        if ts <= count(tsLst) then
            if search(txt[i],"<UserTimeStampString") <> -1 then           
                tagStart=1
                foundUTSS=1
                modTxt=modTxt & txt[i] & newline
            else if search(txt[i],"</UserTimeStampString>") <> -1 then
                tagStart=0                
                modTxt=modTxt & tsLst[ts] & newline 
                ts=ts+1        
            end if
        end if
        
        if tagStart=0 then
            modTxt=modTxt & txt[i] & newline
        end if
           
    end for
  
    if foundUTSS =1 then
        write(vdfPath,modTxt)
    else
        log_error("Could not found tag <UserTimeStampString> in vdf file")
    end if
  	
end createVdf

function toGrib(g:fieldset,ts:number) 
	
	#Save grib file
    grbFile=WORK_DIR & "/tmp_" & ts & ".grib"
    write(grbFile,g)

end toGrib    
    

#====================================================
# Write one timestep of one parameter to vdf/vdc
#====================================================

function toVdf(g:fieldset,varName:string,ts:number,varType:string) 
	
	#Save grib file
    #grbFile=WORK_DIR & "/" & varType & "_" & varName & "_" & ts & ".grib"
    #grbFile=WORK_DIR & "/tmp_" & ts & ".grib"
    grbFile=WORK_DIR & "/tmp.grib"
    write(grbFile,g)
    
    #log_info(" grib: " & varName & " " & count(g))
    
    mv_bin=getenv("METVIEW_BIN")
    
    #Generate the raw binary file
    #rawFile=WORK_DIR & "/" & varType & "_" & varName & "_" & ts & ".raw"
    #rawFile=WORK_DIR & "/tmp_" & ts & ".raw"
    rawFile=WORK_DIR & "/tmp.raw"
    cmd=TORAW_EXE & " " & grbFile & " " & rawFile
    
    ret= shell(cmd)
    if ret <> 0 then
  	   log_error("Conversion of GRIB to raw binary failed!")
  	   log_fail("command: " & cmd)
    end if		
   
	#Add raw binary to vdf/vdc
    cmd=VAPOR_SETUP_CMD & ";" &
        " cd " & OUTPUT_DIR & ";" & 
        " raw2vdf -ts " & ts &
                " -varname " & varName &
                " " & VDF_FILE &
                " " & rawFile
 
    ret= shell(cmd)
    if ret <> 0 then
  	   log_error("Adding raw binary to VAPOR failed!")
  	   log_fail("command: " &cmd)
    end if	                 

end toVdf  
     
function processVar(g:fieldset,varName:string,varType:string)     
       
  for i=1 to count(DATE_LST) do
    
    vdate=DATE_LST[i]["date"]
    vtime=DATE_LST[i]["time"]	
    
    log_info(tab & "--> processing: " & vdate & " " & vtime)
     	
    gRes=nil 	
    loop f in g       
      if grib_get_long(f,"validityDate") = vdate and
         grib_get_long(f,"validityTime") = vtime then
         gRes=gRes & f
      end if
    end loop        
    
    #Write to vdf
    toVdf(gRes,varName,i-1,varType)
    
   end for 
   
end processVar   
 
function prepareElevation()

	log_info("Checking vertical parameters --->")
	
	g=read(data : INPUT_GRB, 
	       param: ELEVATION_PARAM,
	       levtype: LEV_TYPE)
	
	if g = nil then
		#log_fail("Elevation parameter=" &  ELEVATION_PARAM & "not found in input data!")
		if LEV_TYPE  = "ml" and ELEVATION_PARAM = "z" then
		    log_info(tab & "Geopotential is not present! Need to be computed or model levels")
		    computeElevation_ml()
		else
			log_fail("Elevation parameter=" &  ELEVATION_PARAM & " not found in input data!") 
		end if	   
	end if	      
	    
	if ELEVATION_PARAM not in UPPER_PARAMS then
       UPPER_PARAMS = UPPER_PARAMS & [ELEVATION_PARAM]
    end if   
	       
end prepareElevation 
 
function computeElevation_ml()
	
	 t=read(data : INPUT_GRB, 
	       param: "t",
	       levtype: LEV_TYPE)
	    
	 q=read(data : INPUT_GRB, 
	       param: "q",
	       levtype: LEV_TYPE)
	    
	 lnsp=read(data : INPUT_GRB, 
	       param: "lnsp",
	       levtype: LEV_TYPE,
	       level: 1)
	       
	 zs=read(data : INPUT_GRB, 
	       param: "z",
	       levtype: "sfc")      
	       	 
	 if t=nil then 
	 	log_fail("Parameter=t not found on model levels! Height of model levels cannot be computed!")
	 else if q = nil then 
	    log_fail("Parameter=q not found on model levels! Height of model levels cannot be computed!")	 	 
	 else if lnsp = nil then 
	    log_fail("Parameter=lnsp not found! Height of model levels cannot be computed!") 	 
	 else if zs = nil then 
	    log_fail("Surface geopotential is not found! Height of model levels cannot be computed!")
	 end if
	 	 	 
	 gRes=nil 
	 for i=1 to count(DATE_LST) do
	 	
	 	vdate=DATE_LST[i]["date"]
    	vtime=DATE_LST[i]["time"]
 
        tAct=nil 	
        loop f in t        
             if grib_get_long(f,"validityDate") = vdate and
                grib_get_long(f,"validityTime") = vtime then
                tAct=tAct & f
             end if
        end loop
        
        qAct=nil 	
        loop f in q        
             if grib_get_long(f,"validityDate") = vdate and
                grib_get_long(f,"validityTime") = vtime then
                qAct=qAct & f
             end if
        end loop      
          
        lnspAct=nil 	
        loop f in lnsp        
             if grib_get_long(f,"validityDate") = vdate and
                grib_get_long(f,"validityTime") = vtime then
                lnspAct=lnspAct & f
             end if
        end loop 
        
        zsAct=nil 	
        loop f in zs        
             if grib_get_long(f,"validityDate") = vdate and
                grib_get_long(f,"validityTime") = vtime then
                zsAct=zsAct & f
             end if
        end loop 
                       
 	    gRes=gRes & geopotential_on_ml(tAct,qAct,lnspAct,zsAct)
 
    end for

    INPUT_GRB = INPUT_GRB & gRes
    
end computeElevation_ml

 
function nextArg()
  
  CURRENT_ARG_INDEX=CURRENT_ARG_INDEX+1   
  
  if CURRENT_ARG_INDEX <= count(ARGS) then
  	    return ARGS[CURRENT_ARG_INDEX]
  else
  	    return ""
  end if
  	    
end nextArg  

function hasNextArg()
	if CURRENT_ARG_INDEX < count(ARGS) then 
	  return 1
	else 
      return 0
    end if
end hasNextArg         
	
function log_fail(txt:string)
	
	log_error(txt)
	fail()
	
end log_fail	
	
function log_error(txt:string)
	
	if errH=0 then 
		errH=file(WORK_DIR & "/m_err.txt")
	end if		
	write(errH,txt,newline)

end log_error	

function log_info(txt:string)

	if infoH=0 then 
		infoH=file(WORK_DIR & "/m_log.txt")
	end if
	write(infoH,txt,newline)

end log_info


function geopotential_on_ml (T_fs:fieldset, q_fs:fieldset, lnsp_fs:fieldset, zs_fs:fieldset)

    Rd = 287.06
    g  = 9.80665

    #Do some basic checking on the input data - 
    #only check the first field of each fielset for efficiency
    gridtype_t = grib_get_string(T_fs[1],    'gridType')
    gridtype_q = grib_get_string(q_fs[1],    'gridType')
    gridtype_l = grib_get_string(lnsp_fs[1], 'gridType')
    gridtype_z = grib_get_string(zs_fs[1],   'gridType')

    if (gridtype_t = 'sh' or gridtype_q = 'sh' or gridtype_l = 'sh' or gridtype_z = 'sh') then
        log_fail("geopotential_on_ml: fields must be gridded, not spectral")
    end if

    #How many levels are we computing?
    levelSize = count(T_fs) 
	bottomLevel = grib_get_long(T_fs[1], 'level')
    topLevel = grib_get_long(T_fs[levelSize], 'level')

    if bottomLevel < topLevel then
       log_fail("geopotential_on_ml: fields are not sorted in descending order according to model levels")
    end if 
    
    #print("levels: ",bottomLevel," ",topLevel)

    zs   = values(zs_fs)     # extract the data for this field as a vector
    lnsp = values(lnsp_fs)   # extract the data for this field as a vector

    pv = grib_get_double_array(lnsp_fs, 'pv') # for computing pressures at model levels
    sp = exp(lnsp)                          # surface pressure
    z_h = zs                                # orography

    #Get the coefficients for computing the pressures  
    #e.g. L-91: A=1..92, B=93..184
    pv = grib_get_double_array(lnsp_fs, 'pv')
    totalLevNum=count(pv)/2-1
    A = pv[1, totalLevNum+1]
    B = pv[totalLevNum+2, 2*totalLevNum+2]  
    
    if totalLevNum <> bottomLevel then
       log_fail("geopotential_on_ml: bottom level should be " & totalLeveNum & " instead of " & bottomLevel)
    end if   
      
    #A = pv[1, levelSize+1]                #e.g. L-91: A=1..92, B=93..184
    #B = pv[levelSize+2, levelSize+levelSize+2]

    z_f = z_h        # output vector
    z_out_fs = nil   # output fieldset - the end result

    #Compute the bottom pressure (on half-levels)
    #initialise to the lowest numbered level
    Ph_levplusone = (A[totalLevNum+1] + (B[totalLevNum+1] * sp)) 

    # We want to integrate up into the atmosphere, starting at the ground
    # so we start at the lowest level (highest number) and keep
    # accumulating the height as we go.
    # See the IFS documentation:
    # http://www.ecmwf.int/research/ifsdocs/DYNAMICS/Chap2_Discretization4.html
    # For speed and file I/O, we perform the computations with vectors instead
    # of fieldsets.

    for lev=bottomLevel to topLevel by -1 do

		levIdx=bottomLevel-lev+1 #levelSize-lev+1
		
		print("level: ",lev," ",levIdx)
		
        # we assume that the data are in (numerically) descending level order - check if true

        T_level = grib_get_long(T_fs[levIdx], 'level')
        q_level = grib_get_long(q_fs[levIdx], 'level')

        if (T_level <> lev) then
            log_fail ('T field index ', levIdx, ' should be level ', lev, ' but it is ', T_level)
        else if (q_level <> lev) then
            log_fail ('q field index ', levIdx, ' should be level ', lev, ' but it is ', q_level)
        end if


        q = values(q_fs[levIdx]) # extract the data for this field as a vector
        T = values(T_fs[levIdx]) # extract the data for this field as a vector

        # quick check for data consistency

        if ((count(q) <> count(T)) or
            (count(q) <> count(sp)) or
            (count(q) <> count(z_h))) then
            errmsg = 'mvl_geopotential_on_ml: T, Q, LNSP and Z, must have the same number of grid points (they have ' &
                      count(T) & ', ' & count(q) & ', ' & count(sp) & ' and ' & count(z_h) & ' respectively).'
            log_fail(errmsg)
        end if

        T = T*(1.+0.609133*q)   # compute moist temperature
        q = 0                   # free the memory for 'q'

        # compute the pressures (on half-levels)

        Ph_lev = (A[lev] + (B[lev] * sp))

        if lev = 1 then
            dlogP = log(Ph_levplusone/0.1) # Ph[lev] is zero, so don't divide by it
            alpha = log(2)
        else
            dlogP = log(Ph_levplusone/Ph_lev)
            dP    = Ph_levplusone-Ph_lev
            alpha = 1 - Ph_lev/dP*dlogP
        end if

        TRd = T*Rd
        T   = 0      # free memory for 'T'

        # z_f is the geopotential of this full level
        # integrate from previous (lower) half-level z_h to the full level

        z_f = z_h + (TRd*alpha)

        # z_h is the geopotential of 'half-levels'
        # integrate z_h to next half level

        z_h = z_h + (TRd*dlogP)

        # store the result (z_f) in a field and add to the output fieldset
        # (add it to the front, not the end, because we are going 'backwards'
        # through the fields)
        
        out_field = set_values(T_fs[levIdx], z_f) 
        out_field = grib_set_string(out_field, ['shortName', "z"])       
        out_field = grib_set_long(out_field,  ['generatingProcessIdentifier',128,
                                              'level', lev])
             
        z_out_fs = z_out_fs &out_field

        Ph_levplusone  = Ph_lev    # store for the next iteration, will be equivalent to (A[lev+1] + (B[lev+1] * sp))

    end for

    return z_out_fs

end geopotential_on_ml

function listToString(lst)

  s=""
  for i=1 to count(lst) do
  	if i = 1 then
  		s=s & lst[i]
  	else
  		s=s & "," & lst[i]
  	end if
  end for
  
  return s
  
end listToString