1 | ! ************************************************************************** |
---|
2 | |
---|
3 | module scripinterp_mod |
---|
4 | |
---|
5 | ! ========================================================================== |
---|
6 | |
---|
7 | use kinds_mod ! defines common data types |
---|
8 | use constants ! defines common constants |
---|
9 | use iounits ! I/O unit manager |
---|
10 | use netcdf |
---|
11 | use netcdf_mod ! netcdf I/O stuff |
---|
12 | use grids ! module containing grid info |
---|
13 | use remap_vars ! module containing remapping info |
---|
14 | use remap_mod ! module containing remapping routines |
---|
15 | use remap_read ! routines for reading remap files |
---|
16 | |
---|
17 | implicit none |
---|
18 | |
---|
19 | real(kind=dbl_kind), dimension(:), allocatable :: & |
---|
20 | grid_out |
---|
21 | |
---|
22 | integer(kind=int_kind) :: & ! netCDF ids for files and arrays |
---|
23 | ncstat, nc_outfile_id, nc_infile_id |
---|
24 | integer (kind=int_kind), dimension(4) :: & |
---|
25 | input_dims, input_dimids, input_count |
---|
26 | real (kind=dbl_kind), dimension(:), allocatable :: & |
---|
27 | scale |
---|
28 | integer (kind=int_kind), dimension(:), allocatable :: & |
---|
29 | nc_xysize_id, nc_gridsize_id, nc_gridsize, & |
---|
30 | nc_variable_id |
---|
31 | integer :: nc_lon_id, nc_lat_id, nc_array_id |
---|
32 | |
---|
33 | character (char_len) :: & |
---|
34 | map_name ! name for mapping from grid1 to grid2 |
---|
35 | character (1), parameter :: & |
---|
36 | separator = '|' |
---|
37 | |
---|
38 | ! - input namelist variables |
---|
39 | |
---|
40 | character (char_len) :: & |
---|
41 | input_file, & ! filename containing input fields |
---|
42 | interp_file, & ! filename containing remap data (map1) |
---|
43 | input_name ! name of variable to grid |
---|
44 | integer (kind=int_kind), dimension(4) :: & |
---|
45 | input_stride, & ! how much of input array to process |
---|
46 | input_start, & ! where to start |
---|
47 | input_stop ! and where to stop |
---|
48 | character (char_len), dimension(4) :: & |
---|
49 | input_vars ! input variables to be copied |
---|
50 | |
---|
51 | ! - output namelist variables |
---|
52 | |
---|
53 | character (char_len) :: & |
---|
54 | output_file, & ! filename for test results |
---|
55 | output_mode, & ! 'create' or 'append' |
---|
56 | output_name, & ! name of new grid variable |
---|
57 | output_lat, & ! as it says |
---|
58 | output_lon, & ! as it says |
---|
59 | output_ydir ! choose to invert output arrays in y direction |
---|
60 | character (char_len), dimension(4) :: & |
---|
61 | output_dims, & ! name of new grid variable |
---|
62 | output_vars ! variables to copy |
---|
63 | character (char_len), dimension(10) :: & |
---|
64 | output_attributes, & ! attributes of stuff in output file |
---|
65 | output_scaling ! scaling factor to apply to input data to get |
---|
66 | ! correct units in output file |
---|
67 | contains |
---|
68 | |
---|
69 | ! ========================================================================== |
---|
70 | |
---|
71 | subroutine process_grid(nm_in) |
---|
72 | |
---|
73 | !----------------------------------------------------------------------- |
---|
74 | ! - dummy variables |
---|
75 | |
---|
76 | character (char_len) :: & |
---|
77 | nm_in ! name of input namelist file |
---|
78 | |
---|
79 | !----------------------------------------------------------------------- |
---|
80 | ! - local variables |
---|
81 | |
---|
82 | integer (kind=int_kind), dimension(4) :: & |
---|
83 | astart, acount, plus_one |
---|
84 | integer (kind=int_kind), dimension(3) :: & |
---|
85 | write_dims |
---|
86 | integer (kind=int_kind) :: & |
---|
87 | i1, i2, jdim, n, nx, ny, nloop, & |
---|
88 | nc_input_id, nc_input_rank, & |
---|
89 | vstart, vstride, numv |
---|
90 | |
---|
91 | real (kind=dbl_kind), dimension(:), allocatable :: & |
---|
92 | grid1_array |
---|
93 | real (kind=dbl_kind), dimension(:,:), allocatable :: & |
---|
94 | var_out |
---|
95 | |
---|
96 | plus_one(:) = 1 |
---|
97 | |
---|
98 | !----------------------------------------------------------------------- |
---|
99 | |
---|
100 | call read_mappings(nm_in) |
---|
101 | |
---|
102 | !----------------------------------------------------------------------- |
---|
103 | ! - read input grid |
---|
104 | ! - WARNING - lots of assumptions here at the moment |
---|
105 | |
---|
106 | ncstat = nf90_open( input_file, NF90_NOWRITE, nc_infile_id ) |
---|
107 | call netcdf_error_handler(ncstat,"open") |
---|
108 | |
---|
109 | ncstat = nf90_inq_varid( nc_infile_id, input_name, nc_input_id ) |
---|
110 | call netcdf_error_handler(ncstat,"inq_varid") |
---|
111 | |
---|
112 | input_dims(:) = 0 |
---|
113 | ncstat = nf90_inquire_variable( nc_infile_id, nc_input_id, & |
---|
114 | ndims=nc_input_rank, dimids=input_dimids(:) ) |
---|
115 | call netcdf_error_handler(ncstat,"inquire_variable") |
---|
116 | |
---|
117 | do jdim = 1,nc_input_rank |
---|
118 | ncstat = nf90_inquire_dimension(nc_infile_id, & |
---|
119 | input_dimids(jdim), len=input_dims(jdim) ) |
---|
120 | call netcdf_error_handler(ncstat,"inquire_dimension") |
---|
121 | enddo |
---|
122 | |
---|
123 | ! - dimids seem to be returned in storage order so the outer dimension of |
---|
124 | ! - the array as described by the netcdf file becomes the first dimension |
---|
125 | ! - as far as f90 is concerned (confused? you will be!) |
---|
126 | |
---|
127 | do jdim = 1,nc_input_rank |
---|
128 | if (input_stop(jdim) == 0) then |
---|
129 | input_stop(jdim) = input_dims(jdim) |
---|
130 | endif |
---|
131 | input_count(jdim) = input_stop(jdim) - input_start(jdim) + 1 |
---|
132 | enddo |
---|
133 | |
---|
134 | ! - rashly we assume x followed by y |
---|
135 | nx = input_dims(1) |
---|
136 | ny = input_dims(2) |
---|
137 | write(*,fmt='("Input grid dimensions are:",2i6)') nx, ny |
---|
138 | if (nx*ny /= grid1_size) then |
---|
139 | write(6,*) "mismatch between input grid and remap data" |
---|
140 | stop |
---|
141 | endif |
---|
142 | |
---|
143 | ! - calculate number of horizontal slices to process |
---|
144 | ! - at the moment this is not very general and will only work with 3 dimensions |
---|
145 | |
---|
146 | acount(1:nc_input_rank) = & |
---|
147 | (input_stop(1:nc_input_rank)-input_start(1:nc_input_rank)+1) / & |
---|
148 | input_stride(1:nc_input_rank) |
---|
149 | nloop = 1 |
---|
150 | do jdim = 1,nc_input_rank |
---|
151 | nloop = nloop*acount(jdim) |
---|
152 | enddo |
---|
153 | nloop = nloop/grid1_size |
---|
154 | write(6,*) "total slices requested: ",nloop |
---|
155 | |
---|
156 | vstart = input_start(nc_input_rank) ! ie extra var has outer dimension |
---|
157 | vstride = input_stride(nc_input_rank) |
---|
158 | |
---|
159 | ! - in general we cant read in the whole array so do it slice by slice |
---|
160 | ! - slow but sure |
---|
161 | |
---|
162 | write(6,*) "allocating input and output grids" |
---|
163 | allocate( grid1_array(grid1_size)) |
---|
164 | allocate( grid_out(grid2_size) ) |
---|
165 | |
---|
166 | numv = 0 |
---|
167 | do n = 1,4 |
---|
168 | if (trim(input_vars(n)) /= '-' .and. & |
---|
169 | trim(output_vars(n)) /= '-') numv = numv + 1 |
---|
170 | enddo |
---|
171 | |
---|
172 | write_dims(1) = grid2_dims(1) |
---|
173 | write_dims(2) = grid2_dims(2) |
---|
174 | write_dims(3) = nloop |
---|
175 | call define_grid(write_dims(1:3) , 2+numv) |
---|
176 | |
---|
177 | astart(:) = input_start(:) |
---|
178 | astart(3) = astart(3) - input_stride(3) |
---|
179 | acount(:) = 1 |
---|
180 | acount(1) = nx |
---|
181 | acount(2) = ny |
---|
182 | |
---|
183 | do n = 1,nloop |
---|
184 | |
---|
185 | write(6,*) "processing slice: ",n |
---|
186 | astart(3) = astart(3) + input_stride(3) |
---|
187 | ncstat = nf90_get_var(nc_infile_id, nc_input_id, grid1_array, & |
---|
188 | start=astart(1:nc_input_rank), & |
---|
189 | count=acount(1:nc_input_rank)) |
---|
190 | call netcdf_error_handler(ncstat,"get_var") |
---|
191 | |
---|
192 | call calculate_grid(grid1_array, grid_out) |
---|
193 | |
---|
194 | call write_grid(grid_out, n, write_dims(1:2) , 2) |
---|
195 | |
---|
196 | enddo |
---|
197 | |
---|
198 | ! --------------------------------------------------------------------- |
---|
199 | ! - now for any extra variables to copy |
---|
200 | |
---|
201 | if (numv > 0) then |
---|
202 | |
---|
203 | write(6,*) "reading ",numv," extra variables" |
---|
204 | allocate( var_out(nloop,numv) ) |
---|
205 | |
---|
206 | do n = 1,numv |
---|
207 | write(6,*) "looking for variable: ",trim(input_vars(n)) |
---|
208 | ncstat = nf90_inq_varid( nc_infile_id, trim(input_vars(n)), nc_input_id ) |
---|
209 | call netcdf_error_handler(ncstat,"inq_varid") |
---|
210 | |
---|
211 | input_dims(:) = 0 |
---|
212 | ncstat = nf90_inquire_variable( nc_infile_id, nc_input_id, & |
---|
213 | ndims=nc_input_rank, dimids=input_dimids(:) ) |
---|
214 | call netcdf_error_handler(ncstat,"inquire_variable") |
---|
215 | |
---|
216 | if (nc_input_rank /= 1) then |
---|
217 | write(6,*) 'sorry, only rank 1 variables can be copied' |
---|
218 | cycle |
---|
219 | endif |
---|
220 | ncstat = nf90_inquire_dimension(nc_infile_id, & |
---|
221 | input_dimids(1), len=input_dims(1) ) |
---|
222 | call netcdf_error_handler(ncstat,"inquire_dimension") |
---|
223 | |
---|
224 | ncstat = nf90_get_var(nc_infile_id, nc_input_id, var_out(1:nloop,n), & |
---|
225 | start=(/ vstart /), stride=(/ vstride /)) |
---|
226 | call netcdf_error_handler(ncstat,"get_var") |
---|
227 | enddo |
---|
228 | |
---|
229 | call write_extra(var_out, numv+2) |
---|
230 | deallocate(var_out) |
---|
231 | |
---|
232 | endif |
---|
233 | |
---|
234 | ncstat = nf90_close(nc_outfile_id) |
---|
235 | call netcdf_error_handler(ncstat,"out close") |
---|
236 | ncstat = nf90_close(nc_infile_id) |
---|
237 | call netcdf_error_handler(ncstat,"in close") |
---|
238 | |
---|
239 | ! --------------------------------------------------------------------- |
---|
240 | |
---|
241 | deallocate( grid1_array, grid_out) |
---|
242 | |
---|
243 | ! --------------------------------------------------------------------- |
---|
244 | |
---|
245 | end subroutine process_grid |
---|
246 | |
---|
247 | ! ========================================================================== |
---|
248 | |
---|
249 | subroutine define_grid(thedims, therank) |
---|
250 | |
---|
251 | !----------------------------------------------------------------------- |
---|
252 | ! - dummy variables |
---|
253 | |
---|
254 | integer (kind=int_kind) :: & |
---|
255 | therank |
---|
256 | integer (kind=int_kind), dimension(therank) :: & |
---|
257 | thedims |
---|
258 | |
---|
259 | !----------------------------------------------------------------------- |
---|
260 | ! - local variables |
---|
261 | |
---|
262 | integer :: & |
---|
263 | k, n, ilon, ilat, icolon, i1, i2, natt, nvar, id, jd, kd, nd |
---|
264 | character (char_len) :: & |
---|
265 | aname, vname, att |
---|
266 | real (kind=dbl_kind) :: s |
---|
267 | |
---|
268 | ! - netcdf variables |
---|
269 | |
---|
270 | integer :: xtype |
---|
271 | |
---|
272 | !----------------------------------------------------------------------- |
---|
273 | ! - define grid size dimensions |
---|
274 | |
---|
275 | allocate(nc_xysize_id(grid2_rank)) |
---|
276 | allocate(nc_gridsize_id(therank)) |
---|
277 | allocate(nc_gridsize(therank)) |
---|
278 | allocate(nc_variable_id(therank-2)) |
---|
279 | |
---|
280 | !----------------------------------------------------------------------- |
---|
281 | ! - setup a NetCDF file for output |
---|
282 | |
---|
283 | xtype = NF90_FLOAT |
---|
284 | |
---|
285 | write(6,*) 'creating output file' |
---|
286 | ncstat = nf90_create (output_file, NF90_CLOBBER, nc_outfile_id) |
---|
287 | call netcdf_error_handler(ncstat,"create") |
---|
288 | |
---|
289 | write(6,*) 'setting global attributes' |
---|
290 | ncstat = nf90_put_att(nc_outfile_id, NF90_GLOBAL, 'title', map_name) |
---|
291 | call netcdf_error_handler(ncstat,"put_att") |
---|
292 | |
---|
293 | write(6,*) 'setting dimensions' |
---|
294 | do n=1,therank |
---|
295 | if (n .eq. therank .and. therank .gt. 2) then |
---|
296 | write(6,*) ' unlimited dim ',trim(output_dims(n)),' size: ',thedims(n) |
---|
297 | ncstat = nf90_def_dim (nc_outfile_id, output_dims(n), NF90_UNLIMITED, & |
---|
298 | nc_gridsize_id(n)) |
---|
299 | else |
---|
300 | write(6,*) ' dim ',trim(output_dims(n)),' size: ',thedims(n) |
---|
301 | ncstat = nf90_def_dim (nc_outfile_id, output_dims(n), thedims(n), & |
---|
302 | nc_gridsize_id(n)) |
---|
303 | endif |
---|
304 | call netcdf_error_handler(ncstat,"def_dim") |
---|
305 | end do |
---|
306 | nc_gridsize(:) = thedims(1:therank) |
---|
307 | |
---|
308 | ! - at the moment there is an assumption here that the ordering is (lon,lat) |
---|
309 | |
---|
310 | ilon = 1 |
---|
311 | ilat = 2 |
---|
312 | nc_xysize_id(1) = nc_gridsize_id(ilon) |
---|
313 | nc_xysize_id(2) = nc_gridsize_id(ilat) |
---|
314 | |
---|
315 | ! ---------------------------------------------------------------- |
---|
316 | ! - define grid center longitude array |
---|
317 | |
---|
318 | write(6,*) 'defining longitude variable' |
---|
319 | ncstat = nf90_def_var (nc_outfile_id, output_lon, & |
---|
320 | xtype, nc_xysize_id, & |
---|
321 | nc_lon_id) |
---|
322 | call netcdf_error_handler(ncstat,"def_var") |
---|
323 | |
---|
324 | ncstat = nf90_put_att (nc_outfile_id, nc_lon_id, 'units', 'degrees') |
---|
325 | call netcdf_error_handler(ncstat,"put_att") |
---|
326 | |
---|
327 | ! ---------------------------------------------------------------- |
---|
328 | ! - define grid center latitude array |
---|
329 | |
---|
330 | write(6,*) 'defining latitude variable' |
---|
331 | ncstat = nf90_def_var (nc_outfile_id, output_lat, & |
---|
332 | xtype, nc_xysize_id, & |
---|
333 | nc_lat_id) |
---|
334 | call netcdf_error_handler(ncstat,"def_var") |
---|
335 | |
---|
336 | ncstat = nf90_put_att (nc_outfile_id, nc_lat_id, 'units', 'degrees') |
---|
337 | call netcdf_error_handler(ncstat,"put_att") |
---|
338 | |
---|
339 | ! ---------------------------------------------------------------- |
---|
340 | ! - define copy variables array |
---|
341 | |
---|
342 | write(6,*) 'defining copy variables' |
---|
343 | do n = 3,therank |
---|
344 | ncstat = nf90_def_var (nc_outfile_id, output_vars(n-2), & |
---|
345 | xtype, nc_gridsize_id(n), & |
---|
346 | nc_variable_id(n-2)) |
---|
347 | call netcdf_error_handler(ncstat,"def_var") |
---|
348 | enddo |
---|
349 | |
---|
350 | ! ---------------------------------------------------------------- |
---|
351 | ! - define output array |
---|
352 | |
---|
353 | write(6,*) 'defining grid variable' |
---|
354 | ncstat = nf90_def_var (nc_outfile_id, output_name, & |
---|
355 | xtype, nc_gridsize_id, & |
---|
356 | nc_array_id) |
---|
357 | call netcdf_error_handler(ncstat,"def_var") |
---|
358 | |
---|
359 | ! ---------------------------------------------------------------- |
---|
360 | ! - output attributes has to come after all other definitions |
---|
361 | ! - this code currently a bit murky, needs a rewrite |
---|
362 | |
---|
363 | ncstat = nf90_inquire (nc_outfile_id, nVariables=nvar) |
---|
364 | call netcdf_error_handler(ncstat,"inquire") |
---|
365 | do n = 1,10 |
---|
366 | att = trim(output_attributes(n)) |
---|
367 | natt = len(att) |
---|
368 | if (att /= '-') then |
---|
369 | i1 = index(att,separator) |
---|
370 | aname = att(1:i1-1) |
---|
371 | do k = 1,nvar |
---|
372 | ncstat = nf90_inquire_variable(nc_outfile_id, k, vname) |
---|
373 | call netcdf_error_handler(ncstat,"inquire_variable") |
---|
374 | if (vname == aname) then |
---|
375 | i2 = index(att,separator,.true.) |
---|
376 | ncstat = nf90_put_att (nc_outfile_id, k, & |
---|
377 | att(i1+1:i2-1), att(i2+1:natt)) |
---|
378 | call netcdf_error_handler(ncstat,"put_att") |
---|
379 | exit ! from innermost do |
---|
380 | endif |
---|
381 | enddo |
---|
382 | endif |
---|
383 | enddo |
---|
384 | |
---|
385 | ! output scaling |
---|
386 | |
---|
387 | allocate (scale(nvar)) |
---|
388 | scale(:) = 1.0 |
---|
389 | |
---|
390 | do n = 1,10 |
---|
391 | att = trim(output_scaling(n)) |
---|
392 | natt = len(att) |
---|
393 | if (att /= '-') then |
---|
394 | i1 = index(att,separator) |
---|
395 | aname = att(1:i1-1) |
---|
396 | do k = 1,nvar |
---|
397 | ncstat = nf90_inquire_variable(nc_outfile_id, k, vname) |
---|
398 | call netcdf_error_handler(ncstat,"inquire_variable") |
---|
399 | if (vname == aname) then |
---|
400 | i2 = index(att,separator,.true.) |
---|
401 | read(att(i2+1:natt),*) scale(k) |
---|
402 | call netcdf_error_handler(ncstat,"put_att") |
---|
403 | exit ! from innermost do |
---|
404 | endif |
---|
405 | enddo |
---|
406 | endif |
---|
407 | enddo |
---|
408 | |
---|
409 | ! ---------------------------------------------------------------- |
---|
410 | ! - end definition stage |
---|
411 | |
---|
412 | ncstat = nf90_enddef(nc_outfile_id) |
---|
413 | call netcdf_error_handler(ncstat,"enddef") |
---|
414 | |
---|
415 | end subroutine define_grid |
---|
416 | |
---|
417 | ! ========================================================================== |
---|
418 | |
---|
419 | subroutine write_grid(thegrid, thelevel, thedims, therank) |
---|
420 | |
---|
421 | !----------------------------------------------------------------------- |
---|
422 | ! - dummy variables |
---|
423 | |
---|
424 | integer (kind=int_kind), intent(in) :: & |
---|
425 | therank, thelevel |
---|
426 | real (kind=dbl_kind), dimension(:), intent(in) :: & |
---|
427 | thegrid |
---|
428 | integer (kind=int_kind), dimension(therank) :: & |
---|
429 | thedims |
---|
430 | |
---|
431 | !----------------------------------------------------------------------- |
---|
432 | ! - local variables |
---|
433 | |
---|
434 | integer :: & |
---|
435 | k, n, ilon, ilat, icolon, j1, j2, dj, natt, nvar, id, jd, kd, nd |
---|
436 | character (char_len) :: & |
---|
437 | aname, vname, att |
---|
438 | real (kind=dbl_kind), dimension(:,:,:), allocatable :: & |
---|
439 | data |
---|
440 | real (kind=dbl_kind) :: s |
---|
441 | real (kind=dbl_kind), parameter :: todeg = 57.295779513082323 |
---|
442 | integer (kind=int_kind), dimension(3) :: & |
---|
443 | start |
---|
444 | |
---|
445 | ! - netcdf variables |
---|
446 | |
---|
447 | integer :: xtype |
---|
448 | |
---|
449 | !----------------------------------------------------------------------- |
---|
450 | ! - write results to NetCDF file |
---|
451 | |
---|
452 | allocate (data(thedims(1),thedims(2),1)) |
---|
453 | if (output_ydir .eq. 'invert') then |
---|
454 | j1 = thedims(2) |
---|
455 | j2 = 1 |
---|
456 | dj = -1 |
---|
457 | else |
---|
458 | j1 = 1 |
---|
459 | j2 = thedims(2) |
---|
460 | dj = 1 |
---|
461 | endif |
---|
462 | |
---|
463 | if (thelevel .eq. 1) then |
---|
464 | |
---|
465 | ! - grid center latitude array |
---|
466 | |
---|
467 | write(6,*) 'writing latitude variable' |
---|
468 | s = scale(nc_lat_id) |
---|
469 | nd = 0 |
---|
470 | do jd = j1,j2,dj |
---|
471 | do id =1,thedims(1) |
---|
472 | nd = nd + 1 |
---|
473 | data(id,jd,1) = s*todeg*grid2_center_lat(nd) |
---|
474 | enddo |
---|
475 | enddo |
---|
476 | ncstat = nf90_put_var(nc_outfile_id, nc_lat_id, data(:,:,1)) |
---|
477 | call netcdf_error_handler(ncstat,"put_var") |
---|
478 | |
---|
479 | ! - grid center longitude array |
---|
480 | |
---|
481 | write(6,*) 'writing longitude variable' |
---|
482 | s = scale(nc_lon_id) |
---|
483 | nd = 0 |
---|
484 | do jd = j1,j2,dj |
---|
485 | do id =1,thedims(1) |
---|
486 | nd = nd + 1 |
---|
487 | data(id,jd,1) = s*todeg*grid2_center_lon(nd) |
---|
488 | enddo |
---|
489 | enddo |
---|
490 | ncstat = nf90_put_var(nc_outfile_id, nc_lon_id, data(:,:,1)) |
---|
491 | call netcdf_error_handler(ncstat,"put_var") |
---|
492 | |
---|
493 | endif |
---|
494 | |
---|
495 | !----------------------------------------------------------------------- |
---|
496 | ! - new grid |
---|
497 | |
---|
498 | write(6,*) 'writing grid variable' |
---|
499 | n = therank |
---|
500 | s = scale(nc_array_id) |
---|
501 | nd = 0 |
---|
502 | do jd = j1,j2,dj |
---|
503 | do id =1,thedims(1) |
---|
504 | nd = nd + 1 |
---|
505 | data(id,jd,1) = thegrid(nd) |
---|
506 | enddo |
---|
507 | enddo |
---|
508 | write(6,*) 'scaling data ' |
---|
509 | data(:,:,1) = s*data(:,:,1) |
---|
510 | start(:) = (/ 1, 1, thelevel /) |
---|
511 | ncstat = nf90_put_var(nc_outfile_id, nc_array_id, data, start) |
---|
512 | call netcdf_error_handler(ncstat,"put_var") |
---|
513 | deallocate(data) |
---|
514 | |
---|
515 | end subroutine write_grid |
---|
516 | |
---|
517 | ! ========================================================================== |
---|
518 | |
---|
519 | subroutine write_extra(thevars, therank) |
---|
520 | |
---|
521 | real (kind=dbl_kind), dimension(:,:), intent(in) :: & |
---|
522 | thevars |
---|
523 | real (kind=dbl_kind), dimension(:), allocatable :: & |
---|
524 | thedata |
---|
525 | integer (kind=int_kind), intent(in) :: & |
---|
526 | therank |
---|
527 | real (kind=dbl_kind) :: s |
---|
528 | integer :: n |
---|
529 | |
---|
530 | allocate( thedata(size(thevars,1)) ) |
---|
531 | |
---|
532 | ! - copy variable arrays |
---|
533 | |
---|
534 | write(6,*) 'writing copy variables' |
---|
535 | do n = 3,therank |
---|
536 | s = scale(nc_variable_id(n-2)) |
---|
537 | thedata(:) = s*thevars(:,n-2) |
---|
538 | ncstat = nf90_put_var(nc_outfile_id, nc_variable_id(n-2), thedata) |
---|
539 | call netcdf_error_handler(ncstat,"put_var") |
---|
540 | enddo |
---|
541 | |
---|
542 | deallocate( thedata ) |
---|
543 | |
---|
544 | end subroutine write_extra |
---|
545 | |
---|
546 | ! ========================================================================== |
---|
547 | |
---|
548 | subroutine close_grid() |
---|
549 | |
---|
550 | ! close netCDF file |
---|
551 | |
---|
552 | write(6,*) 'closing file' |
---|
553 | ncstat = nf90_close(nc_outfile_id) |
---|
554 | call netcdf_error_handler(ncstat,"close") |
---|
555 | |
---|
556 | end subroutine close_grid |
---|
557 | |
---|
558 | ! ========================================================================== |
---|
559 | |
---|
560 | subroutine read_mappings(nm_in) |
---|
561 | |
---|
562 | !----------------------------------------------------------------------- |
---|
563 | ! - dummy variables |
---|
564 | |
---|
565 | character (char_len) :: & |
---|
566 | nm_in ! name of input namelist file |
---|
567 | |
---|
568 | !----------------------------------------------------------------------- |
---|
569 | ! - local variables |
---|
570 | |
---|
571 | character (char_len) :: & |
---|
572 | dim_name ! netCDF dimension name |
---|
573 | |
---|
574 | integer (kind=int_kind) :: & |
---|
575 | iunit ! unit number for namelist file |
---|
576 | |
---|
577 | !----------------------------------------------------------------------- |
---|
578 | ! - namelist block |
---|
579 | |
---|
580 | namelist /interp_inputs/ input_file, interp_file, input_name, & |
---|
581 | input_stride, input_start, input_stop, & |
---|
582 | input_vars |
---|
583 | namelist /interp_outputs/ output_dims, output_file, output_mode, output_name, & |
---|
584 | output_lat, output_lon, output_ydir, & |
---|
585 | output_scaling, output_vars, output_attributes |
---|
586 | |
---|
587 | !----------------------------------------------------------------------- |
---|
588 | ! - read namelist for file and mapping info |
---|
589 | |
---|
590 | input_stride(:) = 1 |
---|
591 | input_start(:) = 1 |
---|
592 | input_stop(:) = 0 |
---|
593 | output_scaling(:) = '-' |
---|
594 | input_vars(:) = '-' |
---|
595 | output_lon = '-' |
---|
596 | output_lat = '-' |
---|
597 | output_vars(:) = '-' |
---|
598 | output_ydir = 'none' |
---|
599 | output_attributes(:) = '-' |
---|
600 | |
---|
601 | call get_unit(iunit) |
---|
602 | open(iunit, file=nm_in, status='old', form='formatted') |
---|
603 | read(iunit, nml=interp_inputs) |
---|
604 | read(iunit, nml=interp_outputs) |
---|
605 | call release_unit(iunit) |
---|
606 | write(*,nml=interp_inputs) |
---|
607 | write(*,nml=interp_outputs) |
---|
608 | if (trim(output_mode) == "create") then |
---|
609 | if (trim(output_lon) == '-' .or. trim(output_lat) == '-') then |
---|
610 | write(6,*) 'if creating, need to supply lon and lat names' |
---|
611 | stop |
---|
612 | endif |
---|
613 | endif |
---|
614 | |
---|
615 | !----------------------------------------------------------------------- |
---|
616 | ! - read remapping data |
---|
617 | ! - via the scrip package this sets variables: |
---|
618 | ! grid1_size, grid2_size: sizes of input and output grids |
---|
619 | ! grid1_mask, grid2_mask: masks |
---|
620 | ! grid1_rank, grid2_rank: ranks |
---|
621 | |
---|
622 | call read_remap(map_name, interp_file) |
---|
623 | |
---|
624 | end subroutine read_mappings |
---|
625 | |
---|
626 | ! ========================================================================== |
---|
627 | |
---|
628 | subroutine calculate_grid(grid1_array, grid2_array) |
---|
629 | |
---|
630 | !----------------------------------------------------------------------- |
---|
631 | ! - dummy variables |
---|
632 | |
---|
633 | real (kind=dbl_kind), intent(in), dimension(:) :: & |
---|
634 | grid1_array |
---|
635 | real (kind=dbl_kind), intent(out), dimension(:) :: & |
---|
636 | grid2_array |
---|
637 | |
---|
638 | !----------------------------------------------------------------------- |
---|
639 | ! - local variables |
---|
640 | |
---|
641 | integer (kind=int_kind), dimension(:), allocatable :: & |
---|
642 | grid1_imask, grid2_imask, grid2_count |
---|
643 | |
---|
644 | real (kind=dbl_kind), dimension(:), allocatable :: & |
---|
645 | grid1_tmp, & |
---|
646 | grad1_lat, & |
---|
647 | grad1_lon, & |
---|
648 | grad1_latlon, & |
---|
649 | grad1_lat_zero, & |
---|
650 | grad1_lon_zero, & |
---|
651 | grid2_tmp1, & |
---|
652 | grid2_tmp2 |
---|
653 | |
---|
654 | real (kind=dbl_kind) :: & |
---|
655 | delew, delns ! variables for computing bicub gradients |
---|
656 | |
---|
657 | integer (kind=int_kind) :: & |
---|
658 | i,j,n,imin,imax,idiff, & |
---|
659 | ip1,im1,jp1,jm1,nx,ny, & ! for computing bicub gradients |
---|
660 | in,is,ie,iw,ine,inw,ise,isw |
---|
661 | |
---|
662 | logical, parameter :: lat_gradient = .false. |
---|
663 | |
---|
664 | write(6,*) 'starting' |
---|
665 | |
---|
666 | !----------------------------------------------------------------------- |
---|
667 | ! - allocate arrays |
---|
668 | |
---|
669 | allocate (grid1_tmp (grid1_size), & |
---|
670 | grad1_lat (grid1_size), & |
---|
671 | grad1_lon (grid1_size), & |
---|
672 | grad1_lat_zero (grid1_size), & |
---|
673 | grad1_lon_zero (grid1_size), & |
---|
674 | grid1_imask (grid1_size), & |
---|
675 | grid2_tmp1 (grid2_size), & |
---|
676 | grid2_tmp2 (grid2_size), & |
---|
677 | grid2_imask (grid2_size), & |
---|
678 | grid2_count (grid2_size)) |
---|
679 | |
---|
680 | write(6,*) 'allocated' |
---|
681 | write(6,*) grid1_size,grid2_size |
---|
682 | |
---|
683 | grid1_imask(:) = 1 |
---|
684 | grid2_imask(:) = 1 |
---|
685 | where (grid1_mask) |
---|
686 | grid1_imask = 1 |
---|
687 | elsewhere |
---|
688 | grid1_imask = 0 |
---|
689 | endwhere |
---|
690 | where (grid2_mask) |
---|
691 | grid2_imask = 1 |
---|
692 | elsewhere |
---|
693 | grid2_imask = 0 |
---|
694 | endwhere |
---|
695 | |
---|
696 | write(6,*) 'masked' |
---|
697 | |
---|
698 | grad1_lat_zero = zero |
---|
699 | grad1_lon_zero = zero |
---|
700 | nx = input_dims(1) |
---|
701 | ny = input_dims(2) |
---|
702 | write(6,*) nx,ny |
---|
703 | |
---|
704 | !----------------------------------------------------------------------- |
---|
705 | ! - if bicubic, we need 3 gradients in logical space |
---|
706 | |
---|
707 | if (map_type == map_type_bicubic) then |
---|
708 | |
---|
709 | write(6,*) 'bicubic' |
---|
710 | write(6,*) grid1_size |
---|
711 | |
---|
712 | allocate (grad1_latlon (grid1_size)) |
---|
713 | |
---|
714 | do n=1,grid1_size |
---|
715 | |
---|
716 | grad1_lat(n) = zero |
---|
717 | grad1_lon(n) = zero |
---|
718 | grad1_latlon(n) = zero |
---|
719 | |
---|
720 | ! if (n.ge.8000) write(6,*) 0,grid1_mask(n),nx |
---|
721 | if (grid1_mask(n)) then |
---|
722 | |
---|
723 | delew = half |
---|
724 | delns = half |
---|
725 | |
---|
726 | j = (n-1)/nx + 1 |
---|
727 | i = n - (j-1)*nx |
---|
728 | |
---|
729 | ip1 = i+1 |
---|
730 | im1 = i-1 |
---|
731 | jp1 = j+1 |
---|
732 | jm1 = j-1 |
---|
733 | |
---|
734 | if (ip1 > nx) ip1 = ip1 - nx |
---|
735 | if (im1 < 1 ) im1 = nx |
---|
736 | if (jp1 > ny) then |
---|
737 | jp1 = j |
---|
738 | delns = one |
---|
739 | endif |
---|
740 | if (jm1 < 1 ) then |
---|
741 | jm1 = j |
---|
742 | delns = one |
---|
743 | endif |
---|
744 | |
---|
745 | in = (jp1-1)*nx + i |
---|
746 | is = (jm1-1)*nx + i |
---|
747 | ie = (j -1)*nx + ip1 |
---|
748 | iw = (j -1)*nx + im1 |
---|
749 | |
---|
750 | ine = (jp1-1)*nx + ip1 |
---|
751 | inw = (jp1-1)*nx + im1 |
---|
752 | ise = (jm1-1)*nx + ip1 |
---|
753 | isw = (jm1-1)*nx + im1 |
---|
754 | |
---|
755 | ! - compute i-gradient |
---|
756 | |
---|
757 | if (.not. grid1_mask(ie)) then |
---|
758 | ie = n |
---|
759 | delew = one |
---|
760 | endif |
---|
761 | if (.not. grid1_mask(iw)) then |
---|
762 | iw = n |
---|
763 | delew = one |
---|
764 | endif |
---|
765 | |
---|
766 | grad1_lat(n) = delew*(grid1_array(ie) - grid1_array(iw)) |
---|
767 | ! if (n.ge.8000) write(6,*) 1,grad1_lat(n) |
---|
768 | |
---|
769 | ! - compute j-gradient |
---|
770 | |
---|
771 | if (.not. grid1_mask(in)) then |
---|
772 | in = n |
---|
773 | delns = one |
---|
774 | endif |
---|
775 | if (.not. grid1_mask(is)) then |
---|
776 | is = n |
---|
777 | delns = one |
---|
778 | endif |
---|
779 | |
---|
780 | grad1_lon(n) = delns*(grid1_array(in) - grid1_array(is)) |
---|
781 | ! if (n.ge.8000) write(6,*) 2,grad1_lon(n) |
---|
782 | |
---|
783 | ! - compute ij-gradient |
---|
784 | |
---|
785 | delew = half |
---|
786 | if (jp1 == j .or. jm1 == j) then |
---|
787 | delns = one |
---|
788 | else |
---|
789 | delns = half |
---|
790 | endif |
---|
791 | |
---|
792 | if (.not. grid1_mask(ine)) then |
---|
793 | if (in /= n) then |
---|
794 | ine = in |
---|
795 | delew = one |
---|
796 | else if (ie /= n) then |
---|
797 | ine = ie |
---|
798 | inw = iw |
---|
799 | if (inw == n) delew = one |
---|
800 | delns = one |
---|
801 | else |
---|
802 | ine = n |
---|
803 | inw = iw |
---|
804 | delew = one |
---|
805 | delns = one |
---|
806 | endif |
---|
807 | endif |
---|
808 | |
---|
809 | if (.not. grid1_mask(inw)) then |
---|
810 | if (in /= n) then |
---|
811 | inw = in |
---|
812 | delew = one |
---|
813 | else if (iw /= n) then |
---|
814 | inw = iw |
---|
815 | ine = ie |
---|
816 | if (ie == n) delew = one |
---|
817 | delns = one |
---|
818 | else |
---|
819 | inw = n |
---|
820 | ine = ie |
---|
821 | delew = one |
---|
822 | delns = one |
---|
823 | endif |
---|
824 | endif |
---|
825 | |
---|
826 | grad1_lat_zero(n) = delew*(grid1_array(ine) - grid1_array(inw)) |
---|
827 | ! if (n.ge.8000) write(6,*) 3,grad1_lat_zero(n) |
---|
828 | |
---|
829 | if (.not. grid1_mask(ise)) then |
---|
830 | if (is /= n) then |
---|
831 | ise = is |
---|
832 | delew = one |
---|
833 | else if (ie /= n) then |
---|
834 | ise = ie |
---|
835 | isw = iw |
---|
836 | if (isw == n) delew = one |
---|
837 | delns = one |
---|
838 | else |
---|
839 | ise = n |
---|
840 | isw = iw |
---|
841 | delew = one |
---|
842 | delns = one |
---|
843 | endif |
---|
844 | endif |
---|
845 | |
---|
846 | if (.not. grid1_mask(isw)) then |
---|
847 | if (is /= n) then |
---|
848 | isw = is |
---|
849 | delew = one |
---|
850 | else if (iw /= n) then |
---|
851 | isw = iw |
---|
852 | ise = ie |
---|
853 | if (ie == n) delew = one |
---|
854 | delns = one |
---|
855 | else |
---|
856 | isw = n |
---|
857 | ise = ie |
---|
858 | delew = one |
---|
859 | delns = one |
---|
860 | endif |
---|
861 | endif |
---|
862 | |
---|
863 | grad1_lon_zero(n) = delew*(grid1_array(ise) - grid1_array(isw)) |
---|
864 | grad1_latlon(n) = delns*(grad1_lat_zero(n) - grad1_lon_zero(n)) |
---|
865 | ! if (n.ge.8000) write(6,*) 4,grad1_lon_zero(n),grad1_latlon(n) |
---|
866 | |
---|
867 | endif |
---|
868 | enddo |
---|
869 | |
---|
870 | write(6,*) 'remapping' |
---|
871 | call remap(grid2_array, wts_map1, grid2_add_map1, grid1_add_map1, & |
---|
872 | grid1_array, src_grad1=grad1_lat, & |
---|
873 | src_grad2=grad1_lon, src_grad3=grad1_latlon) |
---|
874 | |
---|
875 | print *,'Third order mapping from grid1 to grid2:' |
---|
876 | print *,'----------------------------------------' |
---|
877 | print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array) |
---|
878 | print *,'Grid2 min,max: ',minval(grid2_array ),maxval(grid2_array ) |
---|
879 | |
---|
880 | ! - Conservation Test |
---|
881 | |
---|
882 | print *,'Conservation:' |
---|
883 | print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac) |
---|
884 | print *,'Grid2 Integral = ',sum(grid2_array *grid2_area*grid2_frac) |
---|
885 | |
---|
886 | !----------------------------------------------------------------------- |
---|
887 | ! - a first-order map from grid1 to grid2 |
---|
888 | |
---|
889 | else if (map_type /= map_type_conserv .AND.map_type /= map_type_bicubic) then |
---|
890 | |
---|
891 | write(6,*) 'bilinear or conservative' |
---|
892 | |
---|
893 | call remap(grid2_array, wts_map1, grid2_add_map1, grid1_add_map1,grid1_array) |
---|
894 | |
---|
895 | print *,'First order mapping from grid1 to grid2:' |
---|
896 | print *,'----------------------------------------' |
---|
897 | print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array) |
---|
898 | print *,'Grid2 min,max: ',minval(grid2_array ),maxval(grid2_array ) |
---|
899 | |
---|
900 | ! - Conservation Test |
---|
901 | |
---|
902 | print *,'Conservation:' |
---|
903 | print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac) |
---|
904 | print *,'Grid2 Integral = ',sum(grid2_array *grid2_area*grid2_frac) |
---|
905 | |
---|
906 | !----------------------------------------------------------------------- |
---|
907 | ! - conservative mappings: |
---|
908 | ! - a second-order map from grid1 to grid2 with only lat grads |
---|
909 | |
---|
910 | else if (map_type == map_type_conserv .AND. lat_gradient) then |
---|
911 | |
---|
912 | call remap(grid2_array, wts_map1, grid2_add_map1, grid1_add_map1, & |
---|
913 | grid1_array, src_grad1=grad1_lat,src_grad2=grad1_lon_zero) |
---|
914 | |
---|
915 | select case (norm_opt) |
---|
916 | case (norm_opt_none) |
---|
917 | grid2_tmp2 = grid2_frac*grid2_area |
---|
918 | where (grid2_tmp2 /= zero) |
---|
919 | grid2_array = grid2_array/grid2_tmp2 |
---|
920 | elsewhere |
---|
921 | grid2_array = zero |
---|
922 | end where |
---|
923 | case (norm_opt_frcarea) |
---|
924 | case (norm_opt_dstarea) |
---|
925 | where (grid2_frac /= zero) |
---|
926 | grid2_array = grid2_array/grid2_frac |
---|
927 | elsewhere |
---|
928 | grid2_array = zero |
---|
929 | end where |
---|
930 | end select |
---|
931 | |
---|
932 | print *,'Second order mapping from grid1 to grid2 (lat only):' |
---|
933 | print *,'----------------------------------------' |
---|
934 | print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array) |
---|
935 | print *,'Grid2 min,max: ',minval(grid2_array ),maxval(grid2_array ) |
---|
936 | |
---|
937 | ! - Conservation Test |
---|
938 | |
---|
939 | print *,'Conservation:' |
---|
940 | print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac) |
---|
941 | print *,'Grid2 Integral = ',sum(grid2_array*grid2_area*grid2_frac) |
---|
942 | |
---|
943 | !----------------------------------------------------------------------- |
---|
944 | ! - conservative mappings: |
---|
945 | ! - a second-order map from grid1 to grid2 both gradients |
---|
946 | |
---|
947 | else if (map_type == map_type_conserv .AND..NOT. lat_gradient) then |
---|
948 | |
---|
949 | call remap(grid2_array,wts_map1,grid2_add_map1,grid1_add_map1, & |
---|
950 | grid1_array, src_grad1=grad1_lat,src_grad2=grad1_lon) |
---|
951 | |
---|
952 | select case (norm_opt) |
---|
953 | case (norm_opt_none) |
---|
954 | grid2_tmp2 = grid2_frac*grid2_area |
---|
955 | where (grid2_tmp2 /= zero) |
---|
956 | grid2_array = grid2_array/grid2_tmp2 |
---|
957 | elsewhere |
---|
958 | grid2_array = zero |
---|
959 | end where |
---|
960 | case (norm_opt_frcarea) |
---|
961 | case (norm_opt_dstarea) |
---|
962 | where (grid2_frac /= zero) |
---|
963 | grid2_array = grid2_array/grid2_frac |
---|
964 | elsewhere |
---|
965 | grid2_array = zero |
---|
966 | end where |
---|
967 | end select |
---|
968 | |
---|
969 | print *,'Second order mapping from grid1 to grid2:' |
---|
970 | print *,'-----------------------------------------' |
---|
971 | print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array) |
---|
972 | print *,'Grid2 min,max: ',minval(grid2_array ),maxval(grid2_array ) |
---|
973 | |
---|
974 | ! - Conservation Test |
---|
975 | |
---|
976 | print *,'Conservation:' |
---|
977 | print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac) |
---|
978 | print *,'Grid2 Integral = ',sum(grid2_array*grid2_area*grid2_frac) |
---|
979 | |
---|
980 | endif |
---|
981 | |
---|
982 | !----------------------------------------------------------------------- |
---|
983 | ! calculate some statistics |
---|
984 | |
---|
985 | grid2_count = zero |
---|
986 | grid2_tmp1 = zero |
---|
987 | grid2_tmp2 = zero |
---|
988 | |
---|
989 | print *,'number of sparse matrix entries ',num_links_map1 |
---|
990 | do n=1,num_links_map1 |
---|
991 | grid2_count(grid2_add_map1(n)) = grid2_count(grid2_add_map1(n)) + 1 |
---|
992 | if (wts_map1(1,n) > one .or. wts_map1(1,n) < zero) then |
---|
993 | grid2_tmp1(grid2_add_map1(n)) = grid2_tmp1(grid2_add_map1(n)) + 1 |
---|
994 | grid2_tmp2(grid2_add_map1(n)) = max(abs(wts_map1(1,n)),grid2_tmp2(grid2_add_map1(n)) ) |
---|
995 | endif |
---|
996 | end do |
---|
997 | |
---|
998 | do n=1,grid2_size |
---|
999 | if (grid2_tmp1(n) > zero) print *,n,grid2_tmp2(n) |
---|
1000 | end do |
---|
1001 | |
---|
1002 | imin = minval(grid2_count, mask=(grid2_count > 0)) |
---|
1003 | imax = maxval(grid2_count) |
---|
1004 | idiff = (imax - imin)/10 + 1 |
---|
1005 | print *,'total number of dest cells ',grid2_size |
---|
1006 | print *,'number of cells participating in remap ',count(grid2_count > zero) |
---|
1007 | print *,'min no of entries/row = ',imin |
---|
1008 | print *,'max no of entries/row = ',imax |
---|
1009 | |
---|
1010 | imax = imin + idiff |
---|
1011 | do n=1,10 |
---|
1012 | print *,'num of rows with entries between ',imin,' - ',imax-1, & |
---|
1013 | count(grid2_count >= imin .and. grid2_count < imax) |
---|
1014 | imin = imin + idiff |
---|
1015 | imax = imax + idiff |
---|
1016 | end do |
---|
1017 | |
---|
1018 | !----------------------------------------------------------------------- |
---|
1019 | ! - deallocate arrays |
---|
1020 | |
---|
1021 | deallocate (grid1_tmp, grad1_lat, grad1_lon, & |
---|
1022 | grad1_lat_zero, grad1_lon_zero, grid1_imask, & |
---|
1023 | grid2_tmp1, grid2_tmp2, & |
---|
1024 | grid2_imask, grid2_count) |
---|
1025 | |
---|
1026 | end subroutine calculate_grid |
---|
1027 | |
---|
1028 | ! ========================================================================== |
---|
1029 | |
---|
1030 | end module scripinterp_mod |
---|
1031 | |
---|