Line data Source code
1 : module horizontal_interpolate
2 :
3 : !
4 : ! Modules Used:
5 : !
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use shr_const_mod, only: SHR_CONST_PI
8 : use cam_abortutils, only: endrun
9 : use scamMod, only: single_column
10 : use cam_logfile, only: iulog
11 : implicit none
12 : private
13 : save
14 :
15 : real(r8) :: gw1(1000), gw2(1000)
16 :
17 : public :: xy_interp_init, xy_interp
18 :
19 : contains
20 0 : subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight_distance)
21 : !------------------------------------------------------------------------------------------------------------
22 : ! This program computes weighting functions to map a variable of (im1,jm1) resolution to (im2,jm2) resolution
23 : ! weight_x(im2,im1) is the weighting function for zonal interpolation
24 : ! weight_y(jm2,jm1) is the weighting function for meridional interpolation
25 : !
26 : ! Author: Chih-Chieh (Jack) Chen -- May 2010
27 : !
28 : !------------------------------------------------------------------------------------------------------------
29 : implicit none
30 : integer, intent(in) :: im1, jm1, im2, jm2
31 : logical, intent(in) :: use_flight_distance !.true. = flight distance, .false. = all mixing ratios
32 : real(r8), intent(in) :: lon0(im1), lat0(jm1)
33 : real(r8), intent(out) :: weight_x(im2,im1), weight_y(jm2,jm1)
34 :
35 0 : real(r8) :: lon1(im1), lat1(jm1)
36 0 : real(r8) :: lon2(im2), lat2(jm2)
37 0 : real(r8) :: slon1(im1+1), slon2(im2+1), slat1(jm1+1), slat2(jm2+1)
38 : real(r8) :: x1_west, x1_east, x2_west, x2_east
39 : real(r8) :: y1_south, y1_north, y2_south, y2_north
40 : integer :: i1, j1, i2, j2
41 :
42 0 : weight_x(:,:) = 0.0_r8
43 0 : weight_y(:,:) = 0.0_r8
44 :
45 : ! lon0 & lat0 are longitude & latitude on the source mesh in radians
46 : ! convert lon1, lat1 from radians to degrees
47 0 : lon1(:) = lon0(:)/SHR_CONST_PI*180.0_r8
48 0 : lat1(:) = lat0(:)/SHR_CONST_PI*180.0_r8
49 :
50 : ! set up lon2, lat2 (target mesh), in CAM convention
51 0 : do i2=1,im2
52 0 : lon2(i2) = (float(i2)-1.0_r8)*360.0_r8/float(im2)
53 : enddo
54 0 : do j2=1,jm2
55 0 : lat2(j2) = -90.0_r8+(float(j2)-1.0_r8)*180.0_r8/(float(jm2)-1.0_r8)
56 : enddo
57 :
58 :
59 : ! set up staggered longitudes (cell edges in x)
60 0 : do i1=2,im1
61 0 : slon1(i1) = (lon1(i1-1)+lon1(i1))/2.0_r8
62 : enddo
63 0 : slon1(1) = lon1(1)-(lon1(2)-lon1(1))/2.0_r8
64 0 : slon1(im1+1) = lon1(im1)+(lon1(im1)-lon1(im1-1))/2.0_r8
65 :
66 0 : do i2=2,im2
67 0 : slon2(i2) = (lon2(i2-1)+lon2(i2))/2.0_r8
68 : enddo
69 0 : slon2(1) = lon2(1)-(lon2(2)-lon2(1))/2.0_r8
70 0 : slon2(im2+1) = lon2(im2)+(lon2(im2)-lon2(im2-1))/2.0_r8
71 :
72 : ! set up staggered lattiudes (cell edges in y)
73 0 : slat1(1)=-90.0_r8
74 0 : do j1=2,jm1
75 0 : slat1(j1) = (lat1(j1-1)+lat1(j1))/2.0_r8
76 : enddo
77 0 : slat1(jm1+1)=90.0_r8
78 :
79 0 : slat2(1)=-90.0_r8
80 0 : do j2=2,jm2
81 0 : slat2(j2)=(lat2(j2-1)+lat2(j2))/2.0_r8
82 : enddo
83 0 : slat2(jm2+1)=90.0_r8
84 :
85 : ! compute Guassian weight for two meshes (discrete form of cos(lat).)
86 0 : do j1=1,jm1
87 0 : gw1(j1) = sin(slat1(j1+1)/180.0_r8*SHR_CONST_PI)-sin(slat1(j1)/180.0_r8*SHR_CONST_PI)
88 : enddo
89 :
90 0 : do j2=1,jm2
91 0 : gw2(j2) = sin(slat2(j2+1)/180.0_r8*SHR_CONST_PI)-sin(slat2(j2)/180.0_r8*SHR_CONST_PI)
92 : enddo
93 :
94 :
95 : ! add 360 to slon1 and slon2
96 0 : slon1(:) = slon1(:)+360.0_r8
97 0 : slon2(:) = slon2(:)+360.0_r8
98 :
99 0 : do i2=1,im2
100 :
101 : ! target grid east-west boundaries
102 0 : x2_west=slon2(i2)
103 0 : x2_east=slon2(i2+1)
104 :
105 0 : do i1=1,im1
106 :
107 : ! source grid east-west boundaries
108 0 : x1_west=slon1(i1)
109 0 : x1_east=slon1(i1+1)
110 :
111 : ! check if there is any overlap between the source grid and the target grid
112 : ! if no overlap, then weighting is zero
113 : ! there are three scenarios overlaps can take place
114 0 : if( (x1_west>=x2_west).and.(x1_east<=x2_east) ) then
115 : ! case 1:
116 : ! x1_west x1_east
117 : ! |-------------------|
118 : ! |---------------------------------|
119 : ! x2_west x2_east
120 0 : if(use_flight_distance) then
121 0 : weight_x(i2,i1) = 1.0_r8
122 : else
123 0 : weight_x(i2,i1) = (x1_east-x1_west)/(x2_east-x2_west)
124 : endif
125 0 : elseif ( (x1_west>=x2_west).and.(x1_west<x2_east) ) then
126 : ! case 2:
127 : ! x1_west x1_east
128 : ! |--------------------------------|
129 : ! |---------------------------------|
130 : ! x2_west x2_east
131 0 : if(use_flight_distance) then
132 0 : weight_x(i2,i1) = (x2_east-x1_west)/(x1_east-x1_west)
133 : else
134 0 : weight_x(i2,i1) = (x2_east-x1_west)/(x2_east-x2_west)
135 : endif
136 0 : elseif ( (x1_east>x2_west).and.(x1_east<=x2_east) ) then
137 : ! case 3:
138 : ! x1_west x1_east
139 : ! |--------------------------------|
140 : ! |---------------------------------|
141 : ! x2_west x2_east
142 0 : if(use_flight_distance) then
143 0 : weight_x(i2,i1) = (x1_east-x2_west)/(x1_east-x1_west)
144 : else
145 0 : weight_x(i2,i1) = (x1_east-x2_west)/(x2_east-x2_west)
146 : endif
147 0 : elseif ( (x1_east>x2_east).and.(x1_west<x2_west) ) then
148 : ! case 4:
149 : ! x1_west x1_east
150 : ! |--------------------------------|
151 : ! |--------------------|
152 : ! x2_west x2_east
153 0 : weight_x(i2,i1) = 1._r8
154 : endif
155 :
156 : enddo
157 : enddo
158 :
159 :
160 : ! consider end points
161 0 : if(slon1(im1+1)>slon2(im2+1)) then
162 : ! case 1:
163 : ! slon1(im1) slon1(im1+1) <--- end point
164 : ! |-------------------------|
165 : ! |----------------|......................|
166 : ! slon2(im2) slon2(im2+1) slon2(2) (note: slon2(im2+1) = slon2(1))
167 0 : if(use_flight_distance) then
168 0 : weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon1(im1+1)-slon1(im1))
169 : else
170 0 : weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon2(2)-slon2(1))
171 : endif
172 : endif
173 :
174 0 : if(slon1(im1+1)<slon2(im2+1)) then
175 : ! case 1:
176 : ! slon1(im1) slon1(im1+1) slon1(2) (note: slon1(im1+1) = slon1(1))
177 : ! |-------------------------|.............................|
178 : ! |-------------------------------|
179 : ! slon2(im2) slon2(im2+1) <--- end point
180 0 : if(use_flight_distance) then
181 0 : weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon1(2)-slon1(1))
182 : else
183 0 : weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon2(2)-slon2(1))
184 : endif
185 : endif
186 :
187 :
188 :
189 0 : do j2=1,jm2
190 : ! target grid north-south boundaries
191 0 : y2_south=slat2(j2)
192 0 : y2_north=slat2(j2+1)
193 :
194 0 : do j1=1,jm1
195 :
196 : ! source grid north-south boundaries
197 0 : y1_south=slat1(j1)
198 0 : y1_north=slat1(j1+1)
199 :
200 : ! check if there is any overlap between the source grid and the target grid
201 : ! if no overlap, then weighting is zero
202 : ! there are three scenarios overlaps can take place
203 : ! note: there is Guassian weight to consider in the meridional direction!
204 :
205 0 : if( (y1_south>=y2_south).and.(y1_north<=y2_north) ) then
206 : ! case 1:
207 : ! y1_south y1_north
208 : ! |-------------------|
209 : ! |---------------------------------|
210 : ! y2_south y2_north
211 0 : if(use_flight_distance) then
212 0 : weight_y(j2,j1) = 1.0_r8
213 : else
214 0 : weight_y(j2,j1) = gw1(j1)/gw2(j2)
215 : endif
216 0 : elseif ( (y1_south>=y2_south).and.(y1_south<y2_north) ) then
217 : ! case 2:
218 : ! y1_south y1_north
219 : ! |--------------------------------|
220 : ! |---------------------------------|
221 : ! y2_south y2_north
222 0 : if(use_flight_distance) then
223 0 : weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south)
224 : else
225 0 : weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2)
226 : endif
227 0 : elseif ( (y1_north>y2_south).and.(y1_north<=y2_north) ) then
228 : ! case 3:
229 : ! y1_south y1_north
230 : ! |--------------------------------|
231 : ! |---------------------------------|
232 : ! y2_south y2_north
233 0 : if(use_flight_distance) then
234 0 : weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)
235 : else
236 0 : weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2)
237 : endif
238 0 : elseif ( (y1_north>y2_north).and.(y1_south<y2_south) ) then
239 : ! case 4:
240 : ! y1_south y1_north
241 : ! |--------------------------------|
242 : ! |---------------------|
243 : ! y2_south y2_north
244 0 : if(use_flight_distance) then
245 0 : weight_y(j2,j1) = 1._r8
246 : else
247 0 : weight_y(j2,j1) = gw1(j1)/gw2(j2)
248 : endif
249 : endif
250 :
251 : enddo
252 : enddo
253 :
254 0 : end subroutine xy_interp_init
255 :
256 0 : subroutine xy_interp(im1,jm1,km1,im2,jm2,pcols,ncols,weight_x,weight_y,var_src,var_trg,lons,lats,count_x,count_y,index_x,index_y)
257 : !-------------------------------------------------------------------------------------------------------------
258 : ! This program interpolates var_src(im1,jm1,km1) to var_trg(im2,jm2,km1) based on weighting functions weight_x & weight_y.
259 : !-------------------------------------------------------------------------------------------------------------
260 : implicit none
261 : integer, intent(in) :: im1 ! source number of longitudes
262 : integer, intent(in) :: jm1 ! source number of latitudes
263 : integer, intent(in) :: km1 ! source/target number of levels
264 : integer, intent(in) :: im2 ! target number of longitudes
265 : integer, intent(in) :: jm2 ! target number of latitudes
266 : integer, intent(in) :: pcols
267 : integer, intent(in) :: ncols
268 : real(r8), intent(in) :: weight_x(im2,im1), weight_y(jm2,jm1)
269 : real(r8), intent(in) :: var_src(im1,jm1,km1)
270 : integer, intent(in) :: lons(pcols), lats(pcols)
271 : integer, intent(in) :: count_x(im2), count_y(jm2)
272 : integer, intent(in) :: index_x(im2,im1), index_y(jm2,jm1)
273 : real(r8), intent(out) :: var_trg(pcols,km1)
274 : integer :: n, i1, j1, k1, i2, j2, ii, jj
275 : real(r8) :: sum_x
276 :
277 0 : var_trg(:,:) = 0.0_r8
278 :
279 :
280 0 : do k1=1,km1
281 0 : do n=1,ncols
282 : ! interpolate in x
283 0 : do jj=1,count_y(lats(n))
284 0 : sum_x = 0.0_r8
285 0 : do ii=1,count_x(lons(n))
286 0 : sum_x = sum_x + var_src(index_x(lons(n),ii),index_y(lats(n),jj),k1)* &
287 0 : weight_x(lons(n),index_x(lons(n),ii))
288 : enddo
289 0 : var_trg(n,k1) = var_trg(n,k1)+sum_x*weight_y(lats(n),index_y(lats(n),jj))
290 : enddo
291 : enddo
292 : enddo
293 :
294 0 : end subroutine xy_interp
295 :
296 :
297 : end module horizontal_interpolate
|
