forked from sandrolubis/Sudden-Stratospheric-Warmings
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathssw_events_selection.ncl
118 lines (83 loc) · 2.83 KB
/
ssw_events_selection.ncl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
; Author: Dr. Sandro Lubis (Feb, 2021)
; A simple NCL's code to calculate the central dates of the SSWs
; following Charlton and Polvani (2007).
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
begin
f = addfile("data/U10_zm.nc", "r") ; input data
u = f->u(time|:, {lev|10}, {lat|60}, lon|0)
time = f->time
date = cd_calendar(time, -2)
TIME = cd_calendar(time, -5)
year = TIME(:, 0)
month = TIME(:, 1)
day = TIME(:, 2)
nyear = count_unique_values(year)
year_unique = get_unique_values(year)
ssw_flag = new(dimsizes(time), integer, -999)
ssw_flag = 0
do n = 0, (nyear - 2)
segment_index = ind( (year.eq.year_unique(n) .and. month.ge.11) .or. (year.eq.year_unique(n + 1) .and. month.le.3) )
u_segment = u(segment_index)
do i = 1, (dimsizes(u_segment) - 1)
if (u_segment(i - 1).ge.0 .and. u_segment(i).lt.0) then
date_i = date(segment_index(i))
date_i_index = ind(date.eq.date_i)
ssw_flag(date_i_index) = 1
end if
end do
delete(segment_index)
delete(u_segment)
end do
reversal_index = ind(ssw_flag.eq.1)
do i = 1, (dimsizes(reversal_index) - 1)
delta_index = ind(date.ge.date(reversal_index(i - 1)) .and. date.le.date(reversal_index(i)))
u_delta = u(delta_index)
u_delta_pos = ind(u_delta.ge.0)
if (dimsizes(u_delta_pos).lt.20) then
ssw_flag(reversal_index(i)) = 0
end if
delete(delta_index)
delete(u_delta)
delete(u_delta_pos)
end do
reversal_index_new = ind(ssw_flag.eq.1)
do i = 0, (dimsizes(reversal_index_new) - 1)
if ( month(reversal_index_new(i)).le.3 ) then
year_i = year(reversal_index_new(i))
ssw_to_apr_index = ind( (date.ge.date(reversal_index_new(i))) .and. ( year.eq.year_i .and. month.le.4 ) )
if ( all(u(ssw_to_apr_index).lt.0) ) then
ssw_flag(reversal_index_new(i)) = 0
else
u_i = u(ssw_to_apr_index)
u_i_pos = ind(u_i.ge.0)
if ( dimsizes(u_i_pos).lt.10 ) then
ssw_flag(reversal_index_new(i)) = 0
end if
delete(u_i)
delete(u_i_pos)
end if
delete(ssw_to_apr_index)
end if
if ( month(reversal_index_new(i)).ge.11 ) then
year_i = year(reversal_index_new(i))
ssw_to_apr_index = ind( (date.ge.date(reversal_index_new(i))) .and. ( year.eq.(year_i + 1) .and. month.le.4 ) )
if ( all(u(ssw_to_apr_index).lt.0) ) then
ssw_flag(reversal_index_new(i)) = 0
else
u_i = u(ssw_to_apr_index)
u_i_pos = ind(u_i.ge.0)
if ( dimsizes(u_i_pos).lt.10 ) then
ssw_flag(reversal_index_new(i)) = 0
end if
delete(u_i)
delete(u_i_pos)
end if
delete(ssw_to_apr_index)
end if
end do
ssw_index = ind(ssw_flag.eq.1)
ssw_events = date(ssw_index)
month_string = (/"Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"/)
print(day(ssw_index) + "-" + month_string(month(ssw_index) - 1) + "-" + year(ssw_index))
asciiwrite("ssw_events_erainterim.txt", ssw_events)
end