-
Notifications
You must be signed in to change notification settings - Fork 6
/
PtRainPerc_MergeSA.mv
92 lines (71 loc) · 3.4 KB
/
PtRainPerc_MergeSA.mv
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
# Metview Macro
# **************************** LICENSE START ***********************************
#
# Copyright 2020 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 ************************************
############################
# #
# PtRainCDF_MergeSA.mv #
# #
############################
Perc_list = list() # initialization of the empty variable that will contain all the CDFs in SA fields
for CodeSA = 1 to NumSA do
CodeSASTR = left_pad_number(CodeSA,NumDigSA)
# ------------------------------------------------------- #
# Reading the input files containing the CDFs for each EM #
# ------------------------------------------------------- #
print(" ")
print(" a. Reading SA n. " & CodeSA & " containing the point rainfall CDF for all raw EMs...")
CDF_AllEM_SA = nil
for CodeEM = EnsMemS to 2 do
CodeEMSTR = left_pad_number(CodeEM,NumDigEM)
FileIN = WDIR21 & "/" & BaseDateSTR & BaseTimeSTR & "/" & StepFSTR & "/" & CodeSASTR & "/" & NameFile21 & "_" & CodeSASTR & "_" & CodeEMSTR & ".grib"
CDF_AllEM_SA = CDF_AllEM_SA & read(FileIN)
end for
template = CDF_AllEM_SA[1]
CDF_AllEM_SA = values(CDF_AllEM_SA)
# ------------------------------------- #
# Computing the percentiles for the CDF #
# ------------------------------------- #
print(" b. Computing percentiles...")
Perc_sa = percentile(CDF_AllEM_SA, Perc)
# Defining the list that will contain the global field for the CDFs
Perc_list = Perc_list & list(Perc_sa)
end for
# ------------------------------------------ #
# Merging the SAs containing the percentiles #
# ------------------------------------------ #
print(" ")
print(" c. Merging the " & NumSA & " SAs into a global field")
# Setting the global sample grib
temp = read(FileIN)
BaseDate = number((grib_get(temp[1], ["dataDate"]))[1][1])
BaseTime = number((grib_get(temp[1], ["dataTime"]))[1][1])
Step = number((grib_get(temp[1], ["stepRange"]))[1][1])
FileGL = read(FileGL)
FileGL = grib_set(FileGL, ["class", "od", "stream", "enfo", "type", "pf", "expver", 1, "indicatorOfParameter", int(PP_Code), "level", PP_level, "levtype", PP_levtype, "dataDate", BaseDate, "dataTime", BaseTime, "stepRange", Step])
Perc_AllEM_GL_grib = duplicate(FileGL,NumPerc)
Perc_AllEM_GL_list = values(Perc_AllEM_GL_grib)
# Merging
for CodePERC = 1 to NumPerc do
tempGL_PERC = nil
for CodeSA = 1 to NumSA do
tempSA = Perc_list[CodeSA]
tempSA_PERC = tempSA[CodePERC]
tempGL_PERC = tempGL_PERC & tempSA_PERC
end for
Perc_AllEM_GL_list[CodePERC] = tempGL_PERC
end for
Perc_AllEM_GL_grib = set_values(Perc_AllEM_GL_grib, Perc_AllEM_GL_list)
for i = 1 to NumPerc do
Perc_AllEM_GL_grib[i] = grib_set_long(Perc_AllEM_GL_grib[i], ["paramId", PP_Code, "numberOfForecastsInEnsemble", NumPerc, "perturbationNumber", i])
end for
# Saving the output file
FileOUT = WDIR4 & "/" & BaseDateSTR & BaseTimeSTR & "/" & StepFSTR & "/" & NameFile4 & ".grib"
print(" ")
print("Saving the output grib file in " & FileOUT)
write(FileOUT, Perc_AllEM_GL_grib)