File: mv_flextra_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 (2464 lines) | stat: -rw-r--r-- 63,504 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
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388
2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447
2448
2449
2450
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464
#!/usr/bin/metview
# Metview Macro

# **************************** LICENSE START ***********************************
#
# Copyright 2012 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 ************************************

#---------------------------------------------------------------
# Config for the various datatypes 
#---------------------------------------------------------------

global operConfCore=(
    class: "od",
    hour: [0,1,2,3,4,5,
          6,7,8,9,10,11,
          12,13,14,15,16,17,
          18,19,20,21,22,23],
    type:["an","fc","fc","fc","fc","fc",
          "an","fc","fc","fc","fc","fc",
          "an","an","fc","fc","fc","fc",
          "an","fc","fc","fc","fc","fc"],
    time:[0,0,0,0,0,0,
          6,6,6,6,6,6,
          12,12,12,12,12,12,
          18,18,18,18,18,18],          
    fluxTime:
         [-18,0,0,0,0,0,
          0,6,6,6,6,6,
          6,12,12,12,12,12,
          12,18,18,18,18,18]           
  )

#Delayed cut-off
global operConf_dcda=(
    class: "od",
    hour:  [0,3,6,9,
            12,15,18,21],
    type:  ["an","fc","an","fc",
           "an","fc","an","fc"]  ,      
    time:  [0,0,6,6,
            12,12,18,18],                   
    flux_3:[-18,0,0,6,
            6,12,12,18],
    flux_6:[-18,-18,0,0,
              6,6,12,12],                         
    startDate: 2004-10-01,
    endDate: 2013-11-18,
    fcStream:[0,"oper",6,"dcda",12,"oper",18,"dcda"],
    anStream:[0,"oper",6,"dcda",12,"oper",18,"dcda"]  
)

#Long window 4dvar
global operConf_lwda=(
    class: "od",
    hour:  [0,3,6,9,
            12,15,18,21],
    type:  ["an","fc","an","fc",
            "an","fc","an","fc"],      
    time:  [0,0,6,6,
            12,12,18,18],                   
    flux_3:[-18,0,0,6,
            6,12,12,18],
    flux_6:[-18,-18,0,0,
              6,6,12,12],           
    startDate: 2013-11-10,
    endDate: 2200-01-01,
    fcStream:[0,"oper",6,"lwda",12,"oper",18,"lwda"],
    anStream:[0,"oper",6,"lwda",12,"oper",18,"lwda"]  
)

global eraIntConf(
  class: "ei",
  hour: [0,3,6,9,12,15,18,21],
  type: ["an","fc","fc","fc","an","fc","fc","fc"],
  time: [0,0,0,0,12,12,12,12],
  flux: [-12,0,0,0,0,12,12,12],
  stream: ["oper","oper","oper","oper","oper","oper","oper","oper"]
)  

global confLst=[operConf_dcda,operConf_lwda] #,eraIntConf]

#---------------------------------------
# Define parameters
#---------------------------------------

global INPUT_SOURCE=""
global INPUT_FIELDS=nil

global CURRENT_RET = nil

global INDEX_DATA = nil
global DATE_INDEX = 1
global TIME_INDEX = 2
global STEP_INDEX = 3
global PARAM_INDEX = 4
global LEVTYPE_INDEX = 5
global LEVEL_INDEX = 6

global DO_FILE_REGRID=1
global AREA=""
global GRID=""
global d_OUTPUT=""
global f_AVAILABLE=""
global TOP_LEVEL=1
global TOP_LEVEL_UNITS=""
global TOP_ML=1
global BOTTOM_ML=-1
global DATA_ML_NUM=-1
global FC_EXPVER=1
global AN_EXPVER=1
global DATE_START=0
global TIME_START=0
global DATE_END=0
global TIME_END=0
global DATE=0
global TIME=0
global STEP=0
global surfParams=nil
global surfFluxParams=nil
global surfAnParams=nil
global mlParams=nil
global fluxFilePrefix="sflux_"
global DEBUG=0
global dTMP=creTmpDir()

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

INTERNAL_DEV_TEST = 0 #0 or 1

#Normal mode
if INTERNAL_DEV_TEST = 0 then
    args=arguments()
    
    print("Start to retrieve/prepare data at ",now(),newline,
          "-----------------------------------------------------")
    
    if count(args) = 0 then
        fail("No arguments specified!")
    end if

    printDebug("Arguments: ",args)

    testDef=()
    #In test mode we take the arguments from a test file 
    if type(args[1]) = 'number' then
        if args[1] = 1 then    
            t=readTest(args[2])         
            args=t[1]
            testDef.fluxLst=t[2]
        end if
    end if    
    
    #Run the task
    main(args,testDef)

#Development test mode
else
    #TEST_MODE=1
    t=readTest("/var/tmp/cgr/mv_test1.txt") 
    args=t[1]
    testDef=(fluxLst: t[2])
    main(args,testDef)
    print(args[2])   
end if

print("Finished!")
stop(0)

#=================================================================
#
# The concrete task
#
#=================================================================

function main(args:list,testDef)

    print("pwd=",getenv("PWD"))
    print("dTmp=",dTMP)
   
    #See if it is a test
    doTest=0
    if number(args[1]) = 1 then
        doTest=1
    end if   

    #Get model type
    model=args[2]
    if model <> "flextra" and model <> "flexpart" then 
        fail("Invalid model argument: ", model,"!", newline,
            "Allowed values are: \"flextra\" or \"flexpart\"")
    end if

    #Define parameter lists
    if model = "flextra" then
        mlParams=["u","v","t"]
        surfParams=["sp","z"]
    
    else
        mlParams=["u","v","t","q"]
        surfParams=["sd","msl","tcc","10u","10v","2t","2d","sp",
            "z","lsm","lcc","mcc","hcc","skt",
            "stl1","swvl1"]
        surfFluxParams=["lsp","cp","sshf","ewss","nsss","ssr"]      
        surfAnParams=["sdor","cvl","cvh","sr"]     	
    end if

    #Get main mode
    mode=args[3]
    if mode <> "fc" and mode <> "period" then 
        fail("Invalid first argument: ",mode,"!",newline,
            "Allowed values are: \"fc\" or \"period\"")
    end if

    #Get check mode
    doCheck=0
    if number(args[4]) = 1 then
        doCheck=1
    end if 

 	#Get input source
 	INPUT_SOURCE=args[5]
	if INPUT_SOURCE <> "mars" and INPUT_SOURCE <> "file" then 
        fail("Invalid argument for INPUT_SOURCE: ",INPUT_SOURCE,"!",newline,
            "Allowed values are: \"mars\" or \"file\"")
    end if
    
    #Get outdir and AVAILABLE file path
    d_OUTPUT=args[6]
    f_AVAILABLE=args[7]

    #Get regrid mode
	DO_FILE_REGRID=number(args[8])
	if INPUT_SOURCE = "mars" then
		DO_FILE_REGRID = 0
	end if	

	#Get area and resolution 
	if DO_FILE_REGRID = 1 or INPUT_SOURCE = "mars" then	
    	AREA=parse(args[9],"/")
    	GRID=parse(args[10],"/") 
	end if
	
    #Get top level
    TOP_LEVEL=number(args[11])
    TOP_LEVEL_UNITS=args[12] 
    if TOP_LEVEL_UNITS <> "ML" and TOP_LEVEL_UNITS <> "HPA" then
        fail("Invalid argument for TOP_LEVEL_UNITS: ",TOP_LEVEL_UNITS," !",newline,
            "Allowed values are: \"ML\" or \"HPA\"")
    end if

    deltaStep=0
    expectedNum=0
     
    #Get "fc" mode scpecific params
    if mode = "fc" then
        expectedNum=16       
        if count(args) < expectedNum then
            fail("Invalid number of arguments in \"fc\" mode: ",count(args)," (should be > ",expectedNum,")")
        end if
   
        FC_EXPVER=args[13]
        if model = "flexpart" and FC_EXPVER <> "1" then
            fail("Invalid argument for FC_EXPVER: ",FC_EXPVER," !",newline,
                "It must always be 1 for FLEXPART!")
        end if        
        
        DATE=number(args[14])
        
        if DATE <= 0 then
			DATE =  number(now() - hour(-DATE*24),"yyyymmdd")			
		end if
        
        TIME=number(args[15])
        if TIME >=24 then
            TIME=TIME/100
        end if
            
        STEP=parse(args[16],"/")
    
        print("options:")
        print("  INPUT_SOURCE=",INPUT_SOURCE)        
        print("  doCheck=",doCheck)
        print("  d_OUTPUT=",d_OUTPUT) 
        print("  f_AVAILABLE=",f_AVAILABLE) 
        print("  DATE=",DATE) 
        print("  TIME=",TIME) 
        print("  STEP=",STEP) 
        if DO_FILE_REGRID = 1 or INPUT_SOURCE  = "mars" then       
        	print("  AREA=",AREA) 
        	print("  GRID=",GRID) 
        else
        	print("  AREA=","original") 
        	print("  GRID=","original") 
       	end if
       	 	
        #print("EXPVER=",fExpver) 
    
        if count(STEP) < 2 then
            fail("Step count=",count(STEP)," At least two steps must be specified!")
        end if  
        
        deltaStep=STEP[2]-STEP[1] 
       
    end if

    #Get "period" mode scpecific params
    if mode = "period" then 
        expectedNum=19
        if count(args) < expectedNum then
            fail("Invalid number of arguments in \"period\" mode: ",count(args)," (should be < ",expectedNum,")")
        end if

        FC_EXPVER=args[13]
        AN_EXPVER=args[14]
        
        if model = "flexpart" and FC_EXPVER <> "1" then
            fail("Invalid arguments for FC_EXPVER: ",FC_EXPVER," !",newline,
                "It must always be 1 for FLEXPART!")
        end if    
        if model = "flexpart" and AN_EXPVER <> "1" then
            fail("Invalid arguments for AN_EXPVER: ",AN_EXPVER," !",newline,
                "It must always be 1 for FLEXPART!")
        end if  
                
        DATE_START=number(args[15])
        TIME_START=number(args[16])
        DATE_END=number(args[17])
        TIME_END=number(args[18])
        STEP=number(args[19])  #parse(arg[17],"/")
        deltaStep=STEP       
        
        print("options:")
        print("  doCheck=",doCheck)
        print("  d_OUTPUT=",d_OUTPUT) 
        print("  f_AVAILABLE=",f_AVAILABLE) 
        print("  DATE_START=",DATE_START) 
        print("  TIME_START=",TIME_START)
        print("  DATE_END=",DATE_END) 
        print("  TIME_END=",TIME_END)
        print("  STEP=",STEP) 
        if DO_FILE_REGRID = 1 or INPUT_SOURCE = "mars" then       
        	print("  AREA=",AREA) 
        	print("  GRID=",GRID) 
        else
        	print("  AREA=","original") 
        	print("  GRID=","original") 
       	end if             
    end if    
   
    #---------------------------------------
    # Get input file list
    #---------------------------------------
  
    if INPUT_SOURCE = "file" then
    	for i=expectedNum+1 to count(args) do
    		f=args[i]
    		printDebug(" input file:" & f)
    		if not(exist(f)) then
        		fail("Input file does not exist!",newline,
            		"file: ",f)	
            end if		
    		INPUT_FIELDS=INPUT_FIELDS & read(f)
    	end for
    	
    	if INPUT_FIELDS = nil then
    		fail("Could not read any input fields!")
    	end if	
    	
    	generate_index()
    	
    	if doCheck = 1 and DO_FILE_REGRID = 0 then   		
    		determineInputGrid()    	
    	end if
    	
    end if
    
    #---------------------------------------
    # Check directories
    #---------------------------------------

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

    #------------------------------------------------
    # Generate description lists for the retrievals
    #------------------------------------------------

    descLst=nil
    fluxDescLst=nil
    outLst=nil

    #Forecast mode
    if mode = "fc" then       
        descLst=getFcDesc()
        if surfFluxParams <> nil then
            fluxDescLst=getFcFluxDesc(descLst,deltaStep)
        end if  

    #Period mode
    else if mode = "period" then
        descLst=getPeriodDesc()
        if surfFluxParams <> nil then
            fluxDescLst=getPeriodFluxDesc(descLst,deltaStep)
        end if 
    end if 
 
    if descLst = nil then
        fail("No date/step specified!")
    end if  

    printDebug("descLst:")
    printDebug(" ",descLst)
    printDebug("------------")
    printDebug("fluxDescLst:")
    printDebug(" ",fluxDescLst)
    printDebug("------------")
    
    #Test the flux descriptions
    if doTest = 1 then
        testFluxDesc(testDef.fluxLst,fluxDescLst)
        stop(0)
    end if

    if surfFluxParams <> nil and count(fluxDescLst) <= count(descLst) then
        fail("The number of flux steps (",count(fluxDescLst),
             ") must must be larger than the number of steps(",count(descLst),")!")
    end if  

    #------------------------------------------- 
    # Find out how many model levels we have
    #-------------------------------------------

    getLevelRange(descLst[1])
    print("Model level range:")
    print("  Top model level=", TOP_ML) 
    print("  Bottom model level=", BOTTOM_ML) 
    if BOTTOM_ML = -1 then
        fail("Could not determine number of model levels in input data!")
    end if  

    #---------------------------------
    # Perform the retrievals
    #---------------------------------

    #Check if the data is already in place
    needRet=0
    if doCheck = 1 then 
        print("Check output files") 
        for i=1 to count(descLst) do				
            d=descLst[i]
            print("  file: " & d.fName)
            #If the file is in place we mark it not to be retrieved
            if checkInputFile(d) = 1 then						
                d.ret=0
                descLst[i]=d
                print(tab,"-file is already in place")
            else
                print(tab,"-data needs to be generated")
                needRet=1
            end if  
        end for
    end if    

    if doCheck=0 or needRet=1 then
 
 		print("  file: " & d.fName)
 
     	#In forecast mode we retrive all the data in one go!!!
    	if mode = "fc" and INPUT_SOURCE =  "mars" then 				
			
			print("Retrieve data")			
			retrieveFcData(descLst,fluxDescLst)	
			if CURRENT_RET <> nil then
				fail("Failed to retrieve forecats data from MARS! Retrieval=",CURRENT_RET)
			end if
			if INPUT_FIELDS = nil then
				fail("Failed to retrieve forecats data from MARS!")
			end if		
 			generate_index()			
 		end if
 	
        #First get the fluxes
        if surfFluxParams <> nil then
            print("Read/retrieve flux data:")
            for i=1 to count(fluxDescLst) do
                getFluxForDate(fluxDescLst[i])
            end for
            print("Process flux data:")
            procSurfFlux(fluxDescLst,deltaStep)           
        end if

        print("Read/retrieve and process data:")

        #Then the other fields
        for i=1 to count(descLst) do
        
            d=descLst[i]        
            if d.ret=1 then
            
                #Get the gribs                
                getDataForDate(d)
            	
                if surfFluxParams <> nil then
                    #print(tab,"-append flux to output")
                    fFluxName=dTMP & d.fName & ".flux"
                    if not exist(fFluxName) then
                        fail("Flux result file does not exist: ",fFluxName)
                    end if    
                    g=read(fFluxName)
                    append(dTMP & d.fName,g)
                end if    
            
                #Move output into target directory
                cmd="mv " & dTMP & d.fName & " " & d_OUTPUT & "/"
                print(tab,"-move data to output")
                shell(cmd) 
            end if    
        end for  
    end if

    #--------------------------------------
    # Generate the AVIALABLE file
    #--------------------------------------

    #First check if we need to generate it
    if doCheck = 1 then 
        if exist(f_AVAILABLE) then 
            if checkAvailableFile(descLst) = 1 then		
                print("File AVAILABLE is already in place with the right contents!")
                return
            end if
        end if
    end if			

    #Now we generate it
    print("Generate a new AVAILABLE file!")

    fp=file(f_AVAILABLE)

    write(fp,"DATE      TIME     FILNAME     SPECIFICATIONS",newline)
    write(fp,"YYYYMMDD  HHMMSS",newline)
    write(fp,"_____________________________________________",newline)

    loop it in descLst
        outName=it.fName #the output file name
        outDate=it.date #the validity date of the output file
        #The format must be (i8,1x,i6,2(6x,a18)) in FORTRAN!!
        write(fp,string(outDate,"yyyymmdd")," ",string(outDate,"HHMMSS"),"      ",outName,"      ON DISC",newline)
    end loop

    fp=0
  
end main  
  
#=========================================================
# Generate the read/retrieve  definition for the fc mode.
#=========================================================

function getFcDesc()

    printDebug("getFcDesc -->")
    
    runDate = date(DATE) + hour(TIME)
    fClass="od"
    fType="fc"
    fExpver=FC_EXPVER
    descLst=nil
  
    printDebug("  runDate: ",runDate)
  
    conf=getConf(fClass,date(DATE))
    
    loop ts in STEP
    
        fDate=runDate + hour(ts)
        printDebug(" step: ", ts)
        printDebug(" fc: ", fDate)
    
        d=(
            class:    conf.class,
            type:     fType,
            stream:   fcStreamForRun(conf,TIME), 
            expVer:   fExpver,
            date:     fDate,
            day:      DATE, 
            time:     TIME,
            step:     ts,
            fName:    dataFileName(fDate),
            ret: 1
        )

        printDebug(" d=",d)
        descLst = descLst & [d]

    end loop
  
    printDebug("<-- getFcDesc_mars")
  
    return descLst

end getFcDesc_mars

#=================================================
#
# Main function to generate/check
# FLEXTRA GRIB  files in "period" mode
#
#=================================================

function getPeriodDesc()

    printDebug("getPeriodDesc -->")
    
    fClass="od"
    fExpverFc=FC_EXPVER
    fExpverAn=AN_EXPVER
    #fExpver=1
 
    startDate = date(DATE_START) + hour(TIME_START)
    endDate = date(DATE_END) + hour(TIME_END)
 
    #Find out the latest forecast run date.
    #If we are after 9 am we set the run date to today 00.
    #Otherwise we set the run date to yesterday 12.
    today=now()
    if hour(today) < 9 then
    	today = date(yyyymmdd(today-1)) + hour(12)
    else
        today = date(yyyymmdd(today))	 
    end if 
  
    #We go through the dates from starting date to end date.
    #We count and check the number of dates just for
    #safety!
 
    descLst=nil
    cnt=1

    actDate = startDate
    while actDate <= endDate and cnt < 1000 do
 
        conf=getConf(fClass,actDate)
       
        actHour=hour(actDate)
        printDebug(" date: " & yyyymmdd(actDate) & " time: " & actHour)
        
        idx=find(conf.hour,actHour)
        if idx = nil then
            fail("Retrieving data for hour=",actHour," is not supported")
        end if    
        
        #Past
        if actDate < today then            
                   
            d=(                
                #expVer:   fExpver,
                class:    conf.class,
                type:     conf.type[idx],
                date:     actDate,
                fName:    dataFileName(actDate),
                ret: 1
              )

            #Forecast
            if d.type = "fc" then
                d.expVer=fExpverFc
                d.step=actHour-conf.time[idx]
                runDate=actDate-hour(d.step)
                d.day=yyyymmdd(runDate)  
                d.time=conf.time[idx]
                d.stream=fcStreamForRun(conf,d.time)   
                if conf.time[idx] <> hour(runDate) then
                    fail("Retriving data for run=",runDate," and step=",d.step," is not supported")
                end if  
      
            #Analysis  
            else if d.type = "an" then
                d.expVer=fExpverAn
                d.day=yyyymmdd(actDate)
                d.time=hour(actDate)
                d.step=0
                d.stream=anStreamForTime(conf,d.time) 
            end if 
      
            descLst=descLst & [d]           
            printDebug(" d=",d)
       
        #Future (must be a forecast)
        else   
            d=(                
                expVer:   fExpverFc,
                class:    conf.class,
                type:     "fc",
                date:     actDate,
                fName:    dataFileName(actDate),
                ret: 1
      	    )
      	    runDate=today
            d.step=(actDate-runDate)*24
            #d.step=actHour-conf.time[idx]
            #runDate=actDate-hour(d.step)
            d.day=yyyymmdd(runDate)  
            #d.time=conf.time[idx]
            d.time=hour(runDate)
            d.stream=fcStreamForRun(conf,d.time)   
            
            #if conf.time[idx] <> hour(runDate) then
            #    fail("Retriving data for run=",runDate," and step=",d.step," is not supported")
            #end if
            
            descLst=descLst & [d]           
            printDebug(" d=", d)
        end if 
   
        actDate = actDate  + hour(STEP)
        cnt = cnt +1
      
    end while

    printDebug("<-- getPeriodDesc")
    
    return descLst

end getPeriodDesc


#=========================================================
# Define the FLEXPART flux read/retrieval in forecast mode
#=========================================================

function getFcFluxDesc(inDescLst:list,deltaStep:number)

    printDebug("getFcFluxDesc -->")

    fExpver=1 #in flexpart users cannot control the expver
    fClass="od"
    fType="fc"
    descLst=nil
   
    for i=1 to count(inDescLst) do 
       
        d=inDescLst[i]

        printDebug("  d=",d)    
   
        #The current fc step
        actDate=d.date
        actDay=yyyymmdd(actDate) 
        actHour=hour(actDate)
           
        #Need to be refined
        conf=getConf(fClass,d.date)
       
        actDay=d.day
        actTime=d.time
        actStep=d.step                    
        actRun=date(actDay) + hour(actTime) 
      
        #The previous fc step            
        prevDate=actDate - hour(deltaStep)           
        prevHour=hour(prevDate)             
        prevStep=actStep-deltaStep       
        
        if prevStep < 0 then   
                      
            ft=fluxRunTime(conf,actHour,deltaStep)
                        
            if ft <> actTime then
                actTime=ft
                if actTime < 0 then
                    actTime=-actTime
                    actDay=yyyymmdd(date(actDay)-1)
                    actRun=date(actDay) + hour(actTime)
                end if                                               
             end if                
       
            actRun=date(actDay) + hour(actTime)
            actStep=(actDate - actRun)*24
                      
            #The previous fc step            
            prevDate=actDate - hour(deltaStep)           
            prevHour=hour(prevDate)             
            prevStep=actStep-deltaStep
            
            #If we work from an input file for step 0 we check if a previous forecast is available. If it is 
            #not available we just get the step=0 fields for the fluxes (i.e. the fluxes for step=0 
            #will be 0)
            if INPUT_SOURCE = "file" and d.step = 0 then
            	if indexed_field_find((date: actDay,time: actTime, step: prevStep),1) = nil then
            		
            		actDate=d.date      
            		actDay=yyyymmdd(actDate)
            		actTime=d.time
            		actHour=hour(actDate)         
            		actStep=0
            		
            		prevDate=d.date                  		         
            		prevStep=0
            		
            	end if
           end if 		            
        end if
        
        fStream=fcStreamForRun(conf,actTime) 
           
        if prevStep < 0 then
            fail("Invalid previous forecast step=",prevStep," would be needed for flux computations.")
        end if    
            
        dPrev=(
            stream: fStream,
            class: fClass,
            type: fType,
            expver: fExpver,
            date: prevDate,
            day:  actDay,
            time: actTime,
            step: prevStep,
            fName: fluxFileName(prevDate,1),
            prev: 1
        )
            
        dAct=(
            stream: fStream,
            class: fClass,
            type: fType,
            expver: fExpver,
            date: actDate,
            day:  actDay,
            time: actTime,                
            step: actStep,
            fName: fluxFileName(actDate,0),
            prev: 0
        )
        
        printDebug("  dPrev=",dPrev)
        printDebug("  dAct=",dAct)
        
        #Try to avoid retrieving the very same fields twice
        addPrev=1
        if count(descLst) > 1 then
            dd=descLst[count(descLst)]
            if dd.prev = 0 and sameFluxRetrieval(dd,dPrev) = 1 then
                addPrev=0
            end if
        end if            
                      
        if addPrev = 1 then              
            descLst = descLst & [dPrev]
        end if
        descLst = descLst & [dAct]        
   
    end for

    printDebug("<-- getFcFluxDesc")

    return descLst

end getFcFluxDesc

#=================================================
#
# Define the FLEXPART flux retrievals in period mode
#
#=================================================

function getPeriodFluxDesc(inDescLst:list,deltaStep:number)

    printDebug("getPeriodFluxDesc -->")

    fExpver=1 #in flexpart user cannot control the expver
    fClass="od"
    fType="fc"
    descLst=nil

    #Find out the latest forecast run date.
    #If we are after 9 am we set the run date to today 00.
    #Otherwise we set the run date to yesterday 12.
    today=now()
    if hour(today) < 9 then
    	today = date(yyyymmdd(today-1)) + hour(12)
    else
        today = date(yyyymmdd(today))	 
    end if 
    
    printDebug("  today:",today)

    for i=1 to count(inDescLst) do 
       
        d=inDescLst[i]
   
        printDebug("  ",d)    
   
        #The current fc step
        actDate=d.date
        actDay=yyyymmdd(actDate) 
        actHour=hour(actDate)
           
        #Need to be refined
        conf=getConf(fClass,d.date)
                                  
        if d.type = "an" then                                         
            actDay=yyyymmdd(actDate)           
            actTime=fluxRunTime(conf,actHour,deltaStep)
                       
            if actTime < 0 then
                actTime=-actTime
                actDay=yyyymmdd(date(actDay)-1)                               
            end if                                                  
        
        else if d.type = "fc"  then                   
            actDay=d.day
            actTime=d.time
            actStep=d.step            
           
            ft=fluxRunTime(conf,actHour,deltaStep)
                            
            if ft <> actTime then
                actTime=ft
                if actTime < 0 then
                    actTime=-actTime
                    actDay=yyyymmdd(date(actDay)-1)
                    actRun=date(actDay) + hour(actTime)
                    
                else if d.date >= today and actTime > hour(today) then
                    actTime=d.time
                    if actStep-deltaStep < 0 then
                        actTime=actTime-6
                        if actTime < 0 then
                            actStep=actStep+6
                            actDay=actDay-hour(24)
                            actTime=-actTime
                        else
                            actStep=actStep+6
                        end if
                     end if                                   
                end if 
                                              
            end if 
          
        else
            fail("Unknown mars type=",d.type," was determined for class=",fClass," and date=",d.date)
        end if    
        
        actRun=date(actDay) + hour(actTime)
        actStep=(actDate - actRun)*24    
        
        #Get the stream 
        fStream=fcStreamForRun(conf,actTime)
            
        #The previous fc step            
        prevDate=actDate - hour(deltaStep)           
        prevHour=hour(prevDate)             
        prevStep=actStep-deltaStep
            
        if prevStep < 0 then
            fail("Invalid previous forecast step=",prevStep," would be needed for flux computations.")
        end if    
            
        dPrev=(
            stream: fStream,
            class: fClass,
            type: fType,
            expver: fExpver,
            date: prevDate,
            day:  actDay,
            time: actTime,
            step: prevStep,
            fName: fluxFileName(prevDate,1),
            prev: 1
        )
            
        dAct=(
            stream: fStream,
            class: fClass,
            type: fType,
            expver: fExpver,
            date: actDate,
            day:  actDay,
            time: actTime,                
            step: actStep,
            fName: fluxFileName(actDate,0),
            prev: 0
        )
        
        #Try to avoid retrieving the very same fields twice
        addPrev=1
        if count(descLst) > 1 then
            dd=descLst[count(descLst)]
            if dd.prev = 0 and sameFluxRetrieval(dd,dPrev) = 1 then
                addPrev=0
            end if
        end if            
                      
        if addPrev = 1 then              
            descLst = descLst & [dPrev]
        end if
        descLst = descLst & [dAct]        
   
    end for

    printDebug("<-- getPeriodFluxDesc")

    return descLst

end getPeriodFluxDesc


#====================================================================
#
# Read data for FLEXTRA/FLEXPART GRIB file for a given date/step
#
#====================================================================

function retrieveFcData(descLst,fluxDescLst)

    if descLst = nil then
    	return 
    end if	

    mlCnt=BOTTOM_ML-TOP_ML+1

	mSteps=buildMarsSteps(descLst)

    retCore=( 
        stream : descLst[1].stream,
        type :   descLst[1].type,
        expver : descLst[1].expVer,   
        date : DATE,
        time : TIME,
        step : mSteps,              
		grid : GRID,
        area : AREA
	)
	
    #----------------------	
	# ml
	#----------------------
	print(tab,"-retrieve model level data")

	ret=(retCore,
		param: mlParams & ["etadot"],
		levType: "ml",
		levelist : [TOP_ML,"to",BOTTOM_ML]
 		)
	
	#print(ret)
	
	CURRENT_RET=ret 
	mlGrb=retrieve(ret)
	CURRENT_RET=nil
	
    if type(mlGrb) = "nil"  then
        fail("MARS retrieval failed for model level data! Retrieval=",ret)		
    end if
	
	INPUT_FIELDS = INPUT_FIELDS & mlGrb
	
    #----------------------	
	# surface
	#----------------------
	print(tab,"-retrieve surface data")	
	
	ret=(retCore,
		param: surfParams,
		levType: "sfc"			
 		)
 		
 	printDebug(ret)	
 		
    CURRENT_RET=ret
    surfGrb=retrieve(ret)
    CURRENT_RET=nil
    
    if type(surfGrb) = "nil" then
        fail("MARS retrieval failed for surface data! Retrieval=",ret)		
    end if
    	 
    INPUT_FIELDS = INPUT_FIELDS & surfGrb
	
  	#----------------------	
	# surface fluxes
	#----------------------
	
	if fluxDescLst <> nil then

		print(tab,"-retrieve surface flux data")	

		fluxDt=buildMarsDt(fluxDescLst)
		for i=1 to count(fluxDt) do
		
    		ret=( 
        		stream :  fluxDt[i].stream,
        		type   :  fluxDt[i].type,
        		expver :  fluxDt[i].expVer, 
        		date  : fluxDt[i].date,
				time  : fluxDt[i].time,
				step  : fluxDt[i].step,      	   
				grid  : GRID,
        		area  : AREA,
        		param: surfFluxParams,
				levType: "sfc"
			)		
			
			printDebug("ret=",ret)
			
			CURRENT_RET=ret
			fluxGrb=retrieve(ret)
			CURRENT_RET=nil
			
			if type(fluxGrb) = "nil" then
        		fail("MARS retrieval failed for surface fluxes! Retrieval=",ret)		
    		end if
			
			INPUT_FIELDS = INPUT_FIELDS & fluxGrb
  		end for
  	end if	
  
  	#----------------------	
	# surface analysis
	#----------------------
		 
	if surfAnParams <> nil and find(STEP,0) <> nil then
	
		print(tab,"-retrieve surface analysis data")
       
        retAnCore=retCore
        retAnCore.type="an"
        retAnCore.step=0
 	
        ret=(retAnCore,
            param : surfAnParams,
            levtype : "sfc"           
        )    
      
        printDebug(ret)
      
		CURRENT_RET=ret
		surfAnGrb=retrieve(ret)
 		CURRENT_RET=nil
 		 
	 	if type(surfAnGrb) = "nil" then
        	fail("MARS retrieval failed for surface analysis! Retrieval=",ret)		
    	end if
	 
    	INPUT_FIELDS = INPUT_FIELDS & surfAnGrb	
	end if
	
	
end retrieveFcData


function buildMarsSteps(dLst)

	dt=nil
	for i=1 to count(dLst) do	
		if dLst[i].ret = 1 then	
	 		dt = dt & [dLst[i].step]
	 	end if	
	end for
	
	return dt 	
	
end buildMarsSteps

function buildMarsDt(dLst)

	dt=nil
	for i=1 to count(dLst) do	
			
	 	found =0 	
		
		for j=1 to count(dt) do
			
			if (found = 0 and 
			    	dLst[i].day = dt[j].date and dLst[i].time = dt[j].time and dLst[i].stream = dt[j].stream and 
			    	dLst[i].type = dt[j].type and dLst[i].expVer = dt[j].expVer) then			 			   
			  			    
			  		d=dt[j]			  	
			  		if type(d.step) = "number" then
			  			d.step = [d.step,dLst[i].step]
			  		else	
			  			d.step =  d.step & [dLst[i].step]			  	
			  		end if			
			  			dt[j]=d
			  	
	   				found = 1
	   			end if
	   	end for		
			
	
	    if found = 0 then				
				dt = dt & [(date: dLst[i].day, time: dLst[i].time, step: dLst[i].step,
				            stream: dLst[i].stream, type: dLst[i].type, expVer: dLst[i].expVer) ] 
		end if
		
	
	end for
	
	return dt

end extractDt	


function getDataForDate(d)

	if INPUT_FIELDS = nil then
    	getDataForDate_mars(d)
    else
    	getDataForDate_file(d)	
    end if
    
end  getDataForDate   
    
#====================================================================
#
# Read data for FLEXTRA/FLEXPART GRIB file for a given date/step
#
#====================================================================

function getDataForDate_file(d)

    print("  ",d.date)

    mlCnt=BOTTOM_ML-TOP_ML+1

    retCore=( 
        data: INPUT_FIELDS,
        date : d.day,
        time : d.time,
        step : d.step        
        )
	
	if DO_FILE_REGRID = 1 then
		retCore.grid = GRID
        retCore.area = AREA
	end if
	
    #Retrieve etadot
    ret=(retCore,
        param : "etadot",
        levtype : "ml",
        levelist : [TOP_ML,"to",BOTTOM_ML])
    
    print(tab,"-read model level data")
    
    #etadot=read(ret)
    etadot=indexed_filter(ret)
 
    if type(etadot) = "nil" or count(etadot) <> mlCnt then
        fail("Reading failed for etadot! Request=",ret)		
    end if

    checkDateConsitency(etadot[1],d.date)	

    #Retrieve u, v and t
    ret=(retCore,
        param : mlParams,
        levtype : "ml",
        levelist : [TOP_ML,"to",BOTTOM_ML])
    
    #mlData=read(ret)
    mlData=indexed_filter(ret)
   
    if type(mlData) = "nil" or count(mlData) <> mlCnt*count(mlParams) then
        fail("Reading failed for model levels! Request=",ret)		
    end if

    checkDateConsitency(mlData[1],d.date)

    print(tab,"-read surface data")

    #Retrieve surface parameters
    ret=(retCore,
        param : surfParams,
        levtype : "sfc"             
    )    
    
    #surfData=read(ret)
	surfData=indexed_filter(ret)

    if type(surfData) = "nil" then
        fail("Reading failed for surface params! Request=",ret)		
    end if
 
    if count(surfData) < count(surfParams) then
        print(tab,"requested params: ",surfParams)
        print(tab,"retreived params: ",grib_get_string(surfData,"shortName"))
        fail("Reading failed for surface params! Asked for " & count(surfParams) &
        " params, but retrived only " &  count(surfData) & "!")
    end if

    #Remove "ssrdc" field that can be returned by mars by mistake
    #Having this field flexpart might produce wrong output
    if count(surfData) > count(surfParams) then
        surfDataTmp=nil
        loop f in surfData
            if grib_get_string(f,"shortName") <> "ssrdc" then
                surfDataTmp = surfDataTmp & f
            end if
        end loop
        surfData = surfDataTmp
    end if

    checkDateConsitency(surfData[1],d.date)
 
    anData=nil
    if surfAnParams <> nil and d.step =0 then
        retAnCore=retCore
        retAnCore.type="an"
 	
        ret=(retAnCore,
            param : surfAnParams,
            levtype : "sfc"           
        )    
                    
        #anData=read(ret)
        anData=indexed_filter(ret)
        
        if type(anData) = "nil" then
            fail("Reading failed for surface analysis params! Retrieval=",ret)		
        end if

        if count(anData) < count(surfAnParams) then
            print(tab,"requested params: ",surfAnParams)
            print(tab,"retreived params: ",grib_get_string(anData,"shortName"))
            fail("Reading failed for surface params! Asked for " & count(surfAnParams) &
            " params, but retrived only " &  count(anData) & "!")
        end if
 
        #Remove "ssrdc" field that can be returned by mars by mistake
        #Having this field flexpart might produce wrong output
        if count(anData) > count(surfAnParams) then
            anDataTmp=nil
            loop f in anData
                if grib_get_string(f,"shortName") <> "ssrdc" then
                    anDataTmp = anDataTmp & f
                end if
            end loop
            anData = anDataTmp
        end if

    end if		
 
    print(tab,"-compute etadot")
    #Compute the etadot needed by FLEXTRA
    res=mvl_flextra_etadot(surfData,etadot)

    if res = nil then
        fail("Etadot computation failed!")		
    end if

    #The result
    res=mlData & res & surfData
    if anData <> nil then
        res = res & anData
    end if
 
    #Write the resulting fieldset 
    write(dTMP & d.fName,res)
 
    #Clear fieldsets
    mlData=0
    etadot=0
    surfData=0
    res=0
 
end getDataForDate_file

#===================================================
#
# Get FLEXTRA/FLEXPART GRIB file for a given data/step
#
#===================================================

function getDataForDate_mars(d)

    print("  ",d.date)

    mlCnt=BOTTOM_ML-TOP_ML+1

    retCore=( stream : d.stream,
        type : d.type,
        expver : d.expVer,
        date : d.day,
        time : d.time,
        step : d.step,
        grid : GRID,
        area : AREA
        )
	
    #Retrieve etadot
    ret=(retCore,
        param : "etadot",
        levtype : "ml",
        levelist : [TOP_ML,"to",BOTTOM_ML]
        )
    
    print(tab,"-retrieve model level data")
    
    etadot=retrieve(ret)
 
    if type(etadot) = "nil" or count(etadot) <> mlCnt then
        fail("MARS retrieval failed for etadot! Retrieval=",ret)		
    end if

    checkDateConsitency(etadot[1],d.date)	

    #Retrieve u, v and t
    ret=(retCore,
        param : mlParams,
        levtype : "ml",
        levelist : [TOP_ML,"to",BOTTOM_ML]
       )
    
    mlData=retrieve(ret)
   
    if type(mlData) = "nil" or count(mlData) <> mlCnt*count(mlParams) then
        fail("MARS retrieval failed for model levels! Retrieval=",ret)		
    end if

    checkDateConsitency(mlData[1],d.date)

    print(tab,"-retrieve surface data")

    #Retrieve surface parameters
    ret=(retCore,
        param : surfParams,
        levtype : "sfc"    
    )    
    
    surfData=retrieve(ret)

    if type(surfData) = "nil" then
        fail("MARS retrieval failed for surface params! Retrieval=",ret)		
    end if
 
    if count(surfData) < count(surfParams) then
        print(tab,"requested params: ",surfParams)
        print(tab,"retreived params: ",grib_get_string(surfData,"shortName"))
        fail("MARS retrieval failed for surface params! Asked for " & count(surfParams) &
        " params, but retrived only " &  count(surfData) & "!")
    end if

    #Remove "ssrdc" field that can be returned by mars by mistake
    #Having this field flexpart might produce wrong output
    if count(surfData) > count(surfParams) then
        surfDataTmp=nil
        loop f in surfData
            if grib_get_string(f,"shortName") <> "ssrdc" then
                surfDataTmp = surfDataTmp & f
            end if
        end loop
        surfData = surfDataTmp
    end if

    checkDateConsitency(surfData[1],d.date)
 
    anData=nil
    if surfAnParams <> nil and d.step =0 then
        retAnCore=retCore
        retAnCore.type="an"
 	
        ret=(retAnCore,
            param : surfAnParams,
            levtype : "sfc"
        )    
             
        anData=retrieve(ret)
        
        if type(anData) = "nil" then
            fail("MARS retrieval failed for surface analysis params! Retrieval=",ret)		
        end if

        if count(anData) < count(surfAnParams) then
            print(tab,"requested params: ",surfAnParams)
            print(tab,"retreived params: ",grib_get_string(anData,"shortName"))
            fail("MARS retrieval failed for surface params! Asked for " & count(surfAnParams) &
            " params, but retrived only " &  count(anData) & "!")
        end if
 
        #Remove "ssrdc" field that can be returned by mars by mistake
        #Having this field flexpart might produce wrong output
        if count(anData) > count(surfAnParams) then
            anDataTmp=nil
            loop f in anData
                if grib_get_string(f,"shortName") <> "ssrdc" then
                    anDataTmp = anDataTmp & f
                end if
            end loop
            anData = anDataTmp
        end if

    end if		
 
    print(tab,"-compute etadot")
    #Compute the etadot needed by FLEXTRA
    res=mvl_flextra_etadot(surfData,etadot)

    if res = nil then
        fail("Etadot computation failed!")		
    end if

    #The result
    res=mlData & res & surfData
    if anData <> nil then
        res = res & anData
    end if
 
    #Write the resulting fieldset 
    write(dTMP & d.fName,res)
 
    #Clear fieldsets
    mlData=0
    etadot=0
    surfData=0
    res=0
 
end getDataForDate_mars


function getFluxForDate(d)

    if INPUT_FIELDS = nil then
    	getFluxForDate_mars(d)
    else
    	getFluxForDate_file(d)	
    end if
end  getFluxForDate   
    
#===================================================
#
# Get FLEXTRA/FLEXPART GRIB file for a given date/step
#
#===================================================

function getFluxForDate_mars(d)

    #Retrieve (accumulated) surface fluxes
    if surfFluxParams = nil then
        return
    end if
    
    #We do not need data for for step=0 since it is always 0 for
    #accumulated fluxes
    if d.step=0 then
        return
    end if
    
    print("  ",d.date,"[type=",d.type,",run=",d.time,",step=",d.step,"]")
    
    ret=( stream : d.stream,
        class: d.class,
        type : d.type,
        expver : d.expVer,
        date : d.day,
        time : d.time,
        step : d.step,
        grid : GRID,
        area : AREA,
        param : surfFluxParams,
        levtype : "sfc"      
    )

    sFlux=retrieve(ret)
     
    if type(sFlux) = "nil" or count(sFlux) <> count(surfFluxParams) then
        fail("MARS retrieval failed for fluxes! Retrieval=",ret)		
    end if
    
    #Write the fluxes into a tmp file. The final fluxes will be computed 
    #when all the timesteps were processed
    write(dTMP & d.fName,sFlux)
   
end getFluxForDate_mars

#===================================================
#
# Get FLEXTRA/FLEXPART GRIB file for a given date/step
#
#===================================================

function getFluxForDate_file(d)

    #Retrieve (accumulated) surface fluxes
    if surfFluxParams = nil then
        return
    end if
    
    print("  ",d.date,"[run=",d.time,",step=",d.step,"]")
    
    retCore=( 
        data: INPUT_FIELDS       
    )
	
	if DO_FILE_REGRID = 1 then
		retCore.grid = GRID
        retCore.area = AREA
	end if
    
    ret=(retCore,      
        date : d.day,
        time : d.time,
        step : d.step,    
        param : surfFluxParams,
        levtype : "sfc"
    )

    #sFlux=read(ret)
    sFlux=indexed_filter(ret)
     
    if type(sFlux) = "nil" or count(sFlux) <> count(surfFluxParams) then
        fail("Reading failed for fluxes! Request=",ret)		
    end if
    
    #Write the fluxes into a tmp file. The final fluxes will be computed 
    #when all the timesteps were processed
    write(dTMP & d.fName,sFlux)
    
     
    #print(" generate=", d.fName)
    
end getFluxForDate_file

#===================================================
#
# De-accumulate fluxes
#
#===================================================

function procSurfFlux(fluxLst,deltaStep)
 
    printDebug("procSurfFlux -->")
    print("  De-accumulate fluxes")
 
    if count(fluxLst) < 2 then
        fail("Only ",count(fluxLst)," steps are avaliable! At last 2 would be needed.")
    end if    
 
    if fluxLst[1].prev <> 1 then
        fail("The first flux step must be marked as prev=1")
    end if        
    
    #Loop for the flux parameters
    loop p in surfFluxParams
 
        printDebug(" param: " & p)
   
        #De-accumulate the flux
        g=nil
        gPrev=nil
        for i=1 to count(fluxLst) do
            
            if fluxLst[i].prev <> 1 then
                
                fAct=fluxLst[i]
                fPrev=fluxLst[i-1]
                
                outFileName=dataFileName(fAct.date) & ".flux"
       
                printDebug(" act: ", fAct)
                printDebug(" prev: ", fPrev)
                printDebug(" outfile: ", outFileName)
                 
                g=getFluxField(fAct,p) 
               
                if fPrev.step = 0 then
                    res=g
                else         			  
 			        gPrev=getFluxField(fPrev,p)                        
                    if gPrev = nil then
                        fail("Cannot compute flux! Previous field is missing from file=",fPrev)          
                    end if
 			        
 			        res=g-gPrev 			    	    
 			    end if   
 			    
 			    #m -> mm/h
 			    if p="lsp" or p="cp" then
 			        res=res*1000/deltaStep
 			    #J/m^2 ->W/m^2
 			    else if p="sshf" or p="ewss" or p="nsss" or p="ssr" then
 			        res=res/(3600.*deltaStep)
 			    else
 			        fail("Not able to manage flux parameter=",p)
 			    end if    
 			     			     
 			    append(dTMP & outFileName,res) 			    
            
            #get the prev field
            #else if fluxLst[i].step <> 0 then
            #   gPrev=getFluxField(fluxLst[i],p)            
            end if	   
        end for  	
    end loop
 
    printDebug("<-- procSurfFlux")
 
end procSurfFlux

function determineInputGrid()

	if INPUT_FIELDS = nil then
		return 
	end if	

	f=INPUT_FIELDS[1]

	dx=grib_get_long(f,"iDirectionIncrement")
 	dy=grib_get_long(f,"jDirectionIncrement")	
 	south=grib_get_long(f,"latitudeOfLastGridPoint")
 	north=grib_get_long(f,"latitudeOfFirstGridPoint")
 	west=grib_get_long(f,"longitudeOfFirstGridPoint")
 	east=grib_get_long(f,"longitudeOfLastGridPoint") 	
 	gribEdition=grib_get_long(f,"editionNumber")
 
 	# GRIB 1 - millidegrees
 	if gribEdition = 1 then

    	dx=dx/1000
    	dy=dy/1000
    	south=south/1000
    	north=north/1000
    	west=west/1000
    	east=east/1000
 	
 	# GRIB 2 - microdegrees  
 	else if gribEdition = 2  then 
    	dx=dx/1000000
    	dy=dy/1000000
    	south=south/1000000
    	north=north/1000000
    	
    	if west < 0 then 
       		west=west + 360.
    	end if
    	west=west/1000000
 
    	if east < 0 then 
      		east= east + 360.
    	end if
    	east=east/1000000
 	end if 	
    
    GRID=[dx,dy]
    AREA=[south,west,north,east]
 		
end  checkInputGrid	


#===================================================
#
# Check if a given GRIB data file is 
# already in place.
#
#===================================================

function checkInputFile(d)

 #Find out fileName
 fName=d_OUTPUT & "/" & d.fName 	
	
 #Check if file exits
 if not(exist(fName)) then
 	return 0
 end if	

 grb=read(source: fName,
	  levType: "ml")

 #Check if there are any model levels
 if grb = nil then
      print("   No model level were found!")
      return 0
 end if
 
 #Check level num
 levNumGrb=getLevelNum(grb[1])
 if levNumGrb <> DATA_ML_NUM  then
     print("    Level num does not match! Expected: ",DATA_ML_NUM," found: ",levNumGrb)
     return 0
 end if

 #Check experiment version
 #expGrb=grib_get_string(grb[1],"experimentVersionNumber")
 #if compareExpver(dataExpver,expGrb) then
 #    print("    Expver does not match! Expected: ",dataExpver," found: ",expGrb)
 #    return 0
 #end if

 #Check date, time
 dayGrb=grib_get_long(grb[1],"dataDate")
 if d.day <> dayGrb then  
    print("    Date does not match! Expected: ",d.day," found: ",dayGrb)
    return 0
 end if
 
 hourGrb=grib_get_long(grb[1],"dataTime")
 if 100*d.time <> hourGrb then  
    print("    Time does not match! Expected: ",100*d.time," found: ",hourGrb)
    return 0
 end if
 
 #Check step for forecasts
 if d.type = "fc" then
   stepGrb=grib_get_long(grb[1],"stepRange")
   if d.step <> stepGrb then 
      print("    Step does not match! Expected: ",d.step," found: ",stepGrb)
      return 0
   end if
 end if 
 
 #Check grid type 
 gridGrb=grib_get_string(grb[1],"gridType")
 if gridGrb <> "regular_ll"  then
 	print("    Gridtype does not match! Expected: regular_ll found: ",gridGrb)
 	return 0
 end if
 
 #Get metadata from grib file
 dxGrb=grib_get_long(grb[1],"iDirectionIncrement")
 dyGrb=grib_get_long(grb[1],"jDirectionIncrement")
 southGrb=grib_get_long(grb[1],"latitudeOfLastGridPoint")
 northGrb=grib_get_long(grb[1],"latitudeOfFirstGridPoint")
 westGrb=grib_get_long(grb[1],"longitudeOfFirstGridPoint")
 eastGrb=grib_get_long(grb[1],"longitudeOfLastGridPoint") 
 gribEdition=grib_get_long(grb[1],"editionNumber")
 
 dx=GRID[1]
 dy=GRID[2]
 south=AREA[1]
 north=AREA[3]
 west=AREA[2]
 east=AREA[4]

 # GRIB 1 - millidegrees
 if gribEdition = 1 then

    dx=dx*1000
    dy=dy*1000
    south=south*1000
    north=north*1000
    west=west*1000
    east=east*1000
    
    #Check grid resolution
    if dx <> dxGrb then
        print("    dx doest not match! Expected: ",dx," found: ",dxGrb) 
        return 0 
    end if
  
    if dy <> dyGrb then
        print("    dy doest not match! Expected: ",dy," found: ",dyGrb) 
        return 0 
    end if
  
    #Check area
    if south <> southGrb then
        print("    South does not match! Expected: ",south," found: ",southGrb)
        return 0
    end if
  
    if north <> northGrb then
        print("     North does not match! Expected: ",north," found: ",northGrb)
        return 0
    end if
 
    globalMatch=0
    if AREA[2] = -180 and AREA[4] = 180 and 
       westGrb = -180000 and eastGrb = 179000 then
          globalMatch=1
    end if    
    
    if west <> westGrb and globalMatch = 0 then
       print("    West does not match! Expected: ",west," found: ",westGrb)
       return 0
    end if
 
    if east <> eastGrb and globalMatch = 0 then
        print("     East does not match! Expected: ",east," found: ",eastGrb)
       return 0
    end if
    
 # GRIB 2 - microdegrees  
 else if gribEdition = 2  then 
    
    dx=dx*1000000
    dy=dy*1000000
    south=south*1000000
    north=north*1000000
     
    if west < 0 then 
       west=west + 360.
    end if
    west=west*1000000
 
    if east < 0 then 
      east= east + 360.
    end if
    east=east*1000000
 
    #Check grid resolution
    if dx <> dxGrb then
    	print("    dx doest not match! Expected: ",dx," found: ",dxGrb) 
    	return 0 
 	end if
 	
    if dy <> dyGrb then
    	print("    dy doest not match! Expected: ",dy," found: ",dyGrb) 
    	return 0 
 	end if
 	
    #Check area
    if south <> southGrb then
 		print("    South does not match! Expected: ",south," found: ",southGrb)
 		return 0
	end if
	
 	if north <> northGrb then
 		print("     North does not match! Expected: ",north," found: ",northGrb)
 		return 0
 	end if
 
    globalMatch=0
    if AREA[2] = -180 and AREA[4] = 180 and 
       westGrb = 180000000 and eastGrb = 179000000 then
        globalMatch=1
    end if    
    
    if west <> westGrb and globalMatch = 0 then
 	    print("    West does not match! Expected: ",west," found: ",westGrb)
 	    return 0
    end if
    
    if east <> eastGrb and globalMatch = 0 then
 	    print("     East does not match! Expected: ",east," found: ",eastGrb)
 	    return 0
    end if
 else
 		print("     Invalid GRIB edition: ",editionNumber)
 	    return 0
 end if
 
 return 1
 
end checkInputFile

#================================================
#
# Check if the AVILABLE file is
# correctly describing the data we have.
#
#================================================

function checkAvailableFile(outLst)

 if not(exist(f_AVAILABLE)) then
 	return 0
 end if	
 
 avLst=nil
 lines=read(f_AVAILABLE) 
 for i=1 to count(lines) do
     f=parse(lines[i]," ")
     if count(f) >= 4 then
       if substring(f[3],1,2) = "EN" then 	    
          avLst=avLst & f[3]
       end if
     end if 
 end for
 
 outAvLst=nil
 loop f in outLst 
 	outAvLst=outAvLst & f.fName
 end loop
 
 if avLst = outAvLst then
 	return 1
 else
 	return 0
 end if
 
 return 0

end checkAvailableFile		 	
 
function checkDateConsitency(f,outFileDate)   
	
    if count(f) <> 1 then
		fail("Cannot check date consistency!")
	end if	
	
	valDate=grib_get_long(f,"validityDate")
 	valTime=grib_get_string(f,"validityTime") 

 	#Valtime has the HHMM format but does not 
 	#contain leading zeros! So we need to change it!
 	len=length(valTime)
 	if len = 4 or len = 3 then
     	h=number(valTime)/100
 	else 
 	 	h=number(valTime)
 	end if    
 
 	d = date(valDate) + hour(h)
 	if d <> outFileDate then
        fail("Incosistency between dates!",newline,
            "Expected: ",outFileDate," recived: ",d) 
 	end  if    
end checkDateConsitency
   
function getLevelRange(d)

	f=nil
	if INPUT_SOURCE = "mars" then

		#Retrieve one model level for lnsp
 		f=retrieve(
  			stream : d.stream,
  			type : d.type,
  			expver : d.expVer,
  			date : d.day,
  			time : d.time,
  			step : d.step,
  			param : "lnsp",
  			levtype : "ml",
  			levelist : 1,
  			grid : [1.5/1.5]			
  		)

 		if type(f) = "nil"  then			
 			fail("getLevelRange(): MARS retrieval failed for param: lnsp")	
 		end if
	else
		
		printDebug("count=",count(INPUT_FIELDS))
		loop g in INPUT_FIELDS
			if f = nil then			   
				if grib_get_string(g,"typeOfLevel") = "hybrid" then
					f=g
				end if
			end if
		end loop				
		
		
 		if type(f) = "nil"  then
 			fail("getLevelRange(): Failed to read number of model levels from input file. No model levels were found!")	
 		end if			
	end if		
				
 	#Number of levels
  	pv=grib_get_double_array(f,"pv")
  	nLev=count(pv)/2-1

  	DATA_ML_NUM=nLev

  	if TOP_LEVEL_UNITS = "ML" then
    	TOP_ML=TOP_LEVEL
    	BOTTOM_ML=nLev
    	return
   
  	#Find the topmost model level with larger pressure values than the TOP_LEVEL pressure
  	else 

    	ac=vector(nLev+1)
    	bc=vector(nLev+1)

    	for i=1 to nLev+1 do
      		ac[i]=pv[i]
      		bc[i]=pv[nLev+1+i]
    	end for

    	pres_min=exp(minvalue(f));

    	for i=1 to nLev do
      		if ac[i]+bc[i]*pres_min	> TOP_LEVEL*100 then
        		TOP_ML=i
        		BOTTOM_ML=nLev
        		return
      		end if
    	end for
  end if

end getLevelRange


function getLevelNum(f:fieldset)  
    pv=grib_get_double_array(f,"pv")
    return count(pv)/2-1	
end getLevelNum

function compareExpver(v,vGrib)

  n=number(v)
  if n <> 0 then
      if n < 10 then
           v="000" & v
      else if n < 100 then
           v="00" & v
      else if n < 1000 then
           v="000" & v
      end if
  end if

  if v = vGrib then
       return 0
  else 
       return 1
  end if

end compareExpver

function getFluxField(d,p)
    g=nil
    
    fpath=dTMP & d.fName
    
    if not exist(fpath) then
        fail("Flux file does not exist: ",fpath)
    else
        g=read(source: fpath, param: p)
        if g = nil then
            fail("Cannot read parameter " & p & " from flux file: ",fpath)
        end if	
    end if
    
    if count(g) <> 1 then
        fail("Read ",count(g)," flux fields instead of 1",
             newline,"def=",d)
    end if     
    
    return g
end getFluxField             



function dataFileName(fDate)
    return "EN" & string(fDate,"yymmddHH") 
end dataFileName

function fluxFileName(fDate,prev)
    return fluxFilePrefix & string(fDate,"yymmddHH") & "_" & prev 
end fluxFileName

function sameFluxRetrieval(d1,d2)

   return d1.stream = d2.stream and d1.class = d2.class and d1.type = d2.type and
          d1.expVer = d2.expVer and d1.date = d2.date and d1.day = d2.day and
          d1.time = d2.time and d1.step = d2.step

end sameFluxRetrieval

function printDebug(s)
    if DEBUG = 1 then
        print(s)
    end if
end printDebug    

function printDebug(a,b)
    if DEBUG = 1 then
        print(a,b)
    end if
end printDebug 

function creTmpDir()
    
    d=tmpfile() & ".d"
    shell("mkdir -p " & d)
    return d & "/"
    
end creTmpDir

#==============================================================
# Config methods
#==============================================================

function getConf(classId,actDate)

  loop d in confLst
    if d.class = classId then               
        if actDate  >= d.startDate and actDate <= d.endDate then
            return d           
        end if                    
    end if 
  end loop
    
  fail("Mars class=",classId," is not supported for the specified date=",actDate) 
end getConf

function fluxRunTime(conf,hh,ds)
    
    for i=1 to count(conf.hour) do 
        if conf.hour[i] = hh then
            lst=conf["flux_" & ds]
            if type(lst) = "list" then
                return lst[i]
            end if    
        end if
    end for
        
    fail("No flux time is defined for class=",conf.id," and hour=",hh," and step=",ds)
end  runTimeForHour           
 
 
function fcStreamForRun(conf,runTime)
    
    for i=1 to count(conf.fcStream) by 2 do 
        s=conf.fcStream[i]
        if s = runTime then
            return conf.fcStream[i+1]
        end if
    end for
        
    fail("No mars stream is defined for class=",conf.id," and runTime=",runTime)
end fcStreamForRun         

function anStreamForTime(conf,tt)
    
    for i=1 to count(conf.anStream) by 2 do 
        s=conf.anStream[i]        
        if s = tt then
            return conf.anStream[i+1]
        end if
    end for
        
    fail("No mars stream is defined for class=",conf.id," and analysis time=",tt)
end anStreamForTime


#==============================================================
# Test methods
#==============================================================

function testFluxDesc(testLst,descLst)  
    for i=1 to count(testLst) do        
        t=testLst[i]
        d=descLst[i]        
        if t.date <> d.date then
            fail(testFluxFailed(i,"date",t,d))           
        else if t.day <> d.day then
            fail(testFluxFailed(i,"day",t,d))
        else if t.time <> d.time then
            fail(testFluxFailed(i,"time",t,d)) 
        else if t.step <> d.step then
            fail(testFluxFailed(i,"step",t,d))
        else if t.fName <> d.fName then
            fail(testFluxFailed(i,"fName",t,d))
        else if t.prev <> d.prev then
            fail(testFluxFailed(i,"prev",t,d))
        end if   
    end for
    
end testFluxDesc        

function testFluxFailed(i,key,t,d)
    fail("Test failed at item ",i," ",key," does not match",newline,
         "ntest=",t,newline,"desc=",d)
end testFluxFailed

function readTest(fName)
    
    f=read(fName) 
    t=nil
    
    #Args
    t = t & [parse(f[1],",","string")]

    #Flux desc
    desc=nil   
    d=nil
    f[1]=""
        
    loop r in f

        if r = "_," then                       
            if d <> nil then
                desc = desc & [d]
            end if                
            d=()
            
        else if r <> "" then           
            v=parse(r,"=","string")          
            if hasKey(v,"date") = 1 then                
                d.date=date(v[2])
            else if hasKey(v,"day") = 1 then
                d.day=number(v[2])              
            else if hasKey(v,"time") = 1 then
                d.time=number(v[2])   
            else if hasKey(v,"step") = 1 then
                d.step=number(v[2])  
            else if hasKey(v,"fName") = 1 then
                d.fName=v[2]  
            else if hasKey(v,"prev") = 1 then
                d.prev=number(v[2])              
            end if
       end if
     end loop
    
     t=t & [desc]
     return t
      
end readTest             
     
function hasKey(v,key)      
    
    if count(v) <> 2 then
        return 0
    end if  
    
    if search(v[1],key) <> -1 then
        return 1
    end if
    return 0
end hasKey        

function generate_index()

	INDEX_DATA = grib_get(INPUT_FIELDS,["date:l","time:l","step:l","shortName","mars.levtype","level"])

	#print(INDEX_DATA)

end generate_index

# An extrememly simple filter adapted to the needs of 
# this macro!
function indexed_filter(d)

	indexLst=indexed_field_find(d,0)
	f=nil	
	if indexLst <> nil then 
	
		for i=1 to count(indexLst) do
			f= f & INPUT_FIELDS[indexLst[i]]
		end for
		
		if DO_FILE_REGRID = 1 then
    		f=f+0    		
    		f=read(data: f, area: d.area, grid: d.grid)    	
    	end if			
	end if	
	
	return f			
	
end indexed_filter

function indexed_field_find(d,checkExistanceOnly)

	indexLst=nil
	
	param = nil
	if type(d.param) = "list" then
		param = d.param
	else if type(d.param) <> "nil" then
		param = [d.param]		
	end if
	
	levType = nil
	if type(d.levtype) = "list" then
		levType = d.levtype
	else if type(d.levtype) <> "nil" then	
		levType = [d.levtype]	
	end if
	
	level = nil
	if type(d.levelist) = "list" then		
		level= expand_list(d.levelist)
	else if type(d.levelist) <> "nil" then	
		level = [d.levelist]	
	end if
	
	printDebug(" param=",param)
	printDebug(" levtype=",levType)
	printDebug(" level=",level)
	
	printDebug("d=",d)
	
	d.time=d.time*100
	
	for i=1 to count(INPUT_FIELDS) do
		
		#if i < 100 then
		#	print(i,INDEX_DATA[i])
		#end if
		
		if INDEX_DATA[i][DATE_INDEX] = d.date and  
		   INDEX_DATA[i][TIME_INDEX] = d.time and
		   INDEX_DATA[i][STEP_INDEX] = d.step then
				
				b = 1
		   		if param <> nil then
		   			b = (INDEX_DATA[i][PARAM_INDEX] in param )	
		   		end if
		   		
		   		if b = 1 then			   			
		   			if levType <> nil then
		   				b = (INDEX_DATA[i][LEVTYPE_INDEX] in levType) 
		   			end if
		   		end if	
		  
		   		if b = 1 then		   			
		   			if level <> nil then
		   				lval=INDEX_DATA[i][LEVEL_INDEX]
		   				if type(lval) = "string" then
		   					lval=number(lval)
		   				end if	
		   				b = (lval in level)
		   			end if
		   		end if	
		   	
		   		if b = 1 then
		   		  indexLst = indexLst & [i]
		   		  if checkExistanceOnly  = 1 then
		   		  	   return indexLst
		   		  end if	
		   		end if  
		end if
    end for
    
    return indexLst
		   
end indexed_field_find



#expand lists in the format of [NUM1,"TO",NUM2]
function expand_list(lst)

	res=nil
	if count(lst) = 3 then
		if lowercase(lst[2]) = "to" then
			n1=lst[1]
			n2=lst[3]
			if n1 = "string" then
				n1=number(n1)
			end if
			if n2 = "string" then
				n2=number(n2)
			end if
						
			for i=n1 to n2 do
				res = res & [i]
			end for
		end if
	end if
				
	if res = nil then
		res = lst
	end if
	
	return res	

end expand_list