source: trunk/northwind.pro @ 6

Last change on this file since 6 was 2, checked in by pinsard, 18 years ago

initial import from /usr/work/fvi/OPA/geomag/

File size: 3.5 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;
5; NAME: northwind.pro
6;
7; PURPOSE: Check on in the intput field the level of coherence of winds nearby north pole,
8;          and corrects it. Thi stask is performed before interpolation
9;
10; CATEGORY: Subroutine
11;
12; CALLING SEQUENCE: northwind
13;
14; INPUTS:
15;
16;
17; KEYWORD PARAMETERS: None
18;
19; OUTPUTS:
20;          zresul : the same wind field with its upward stripe
21;          corrected in order to get coherence of winds on north polte
22;
23; COMMON BLOCKS:
24;          common_interp.pro
25;
26; SIDE EFFECTS:
27;
28; RESTRICTIONS:
29;
30; EXAMPLE:
31;
32; MODIFICATION HISTORY:
33;
34;------------------------------------------------------------
35;------------------------------------------------------------
36
37
38pro northwind, datglo,zdata_name,t
39
40@common_interp
41
42; Treatment of north pole stripe
43
44   
45; Get from global array, the last stripe and the last but one stripe
46 
47;    Initialising data arrays
48
49     datglo_rect  = replicate(0.,2,jpiatm)
50     datglo_rectb = replicate(0.,2,jpiatm)
51
52; Put data of the last two stripes into arrays
53
54     
55     datglo_rect(0,*) = datglo(*,jpjatm-1,0)
56     datglo_rect(1,*) = datglo(*,jpjatm-1,1)
57
58     datglo_rectb(0,*)= datglo(*,jpjatm-2,0)
59     datglo_rectb(1,*)= datglo(*,jpjatm-2,1)
60
61     
62
63; Converts data of these two lines into polar coordinates
64
65     add_angle=0.
66
67     if (zdata_name EQ data_u_name) then add_angle=!PI/2
68
69
70     datglo_polar=CV_COORD(FROM_RECT=datglo_rect,/TO_POLAR)
71     datglo_polar(0,*)=datglo_polar(0,*)+2*!PI*indgen(jpiatm)/jpiatm
72
73     datglo_polarb=CV_COORD(FROM_RECT=datglo_rect,/TO_POLAR)
74     datglo_polarb(0,*)=datglo_polar(0,*)+2*!PI*indgen(jpiatm)/jpiatm
75   
76
77; Coherence of the last line needs to be corrected
78
79     print, 'CHECKING LAST LINE COHERENCE...'
80
81     mean_mod=replicate(total(datglo_polar(1,*))/jpiatm,jpiatm)
82     mean_arg=replicate((total(datglo_polar(0,*))/jpiatm+2*!PI*west_lon/360.+!PI*(jpiatm-1.)/float(jpiatm)) ,jpiatm)
83     mean_argb=replicate((total(datglo_polarb(0,*))/jpiatm+2*!PI*west_lon/360.+!PI*(jpiatm-1.)/float(jpiatm)) ,jpiatm)
84
85     
86
87     ; Checks whether module is about the same on the North Pole line :
88   
89     test_mod=abs(mean_mod-datglo_polar(1,*))
90
91     IF ( max(test_mod)/mean_mod(0) ) GT 0.01 THEN BEGIN
92   
93          print, 'MODULE IS NOT COHERENT AND WILL BE CORRECTED TO A MEAN VALUE OF :', mean_mod(0)
94         
95          datglo_polar(1,*)=mean_mod
96         
97     ENDIF
98   
99
100     ; Checks whether argument is coherent on the North Pole line :
101
102     test_arg=abs(mean_arg-datglo_polar(0,*))
103
104     IF ( max(test_arg)/mean_arg(0)) GT 0.01 THEN BEGIN
105
106          print, 'ARGUMENT IS NOT COHERENT AND WILL BE CORRECTED TO A MEAN VALUE :', mean_argb(0)
107 
108          datglo_polar(0,*)=mean_argb-!PI/2-(2*!PI*indgen(jpiatm)/float(jpiatm)+replicate(2.*!PI*west_lon/360.,jpiatm))
109           
110     
111
112     ; Cross Corellation Method for wind coherence on North Pole
113
114     all_phases=indgen(jpiatm)
115     
116     Ux=mean_mod*cos(all_phases*2.*!PI/jpiatm)
117     Uy=mean_mod*sin(all_phases*2.*!PI/jpiatm)
118
119     phyx=C_CORRELATE(Ux,datglo(*,jpjatm-2,0),all_phases)
120
121     phyy=C_CORRELATE(Uy,datglo(*,jpjatm-2,1),all_phases)
122
123
124     deph_x=where(phyx EQ max(phyx))
125     deph_y=where(phyy EQ max(phyy))
126     
127     deph_x=deph_x(0)
128     deph_y=deph_y(0)
129
130     
131
132
133     Ux=shift(Ux,deph_x)
134     Uy=shift(Uy,deph_y)
135
136     datglo(*,jpjatm-1,0)=Ux
137     datglo(*,jpjatm-1,1)=Uy
138 
139
140     ENDIF
141
142 
143   
144
145end
Note: See TracBrowser for help on using the repository browser.