source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/SCRIP1.4/source/scrip_test_repeat.f @ 2352

Last change on this file since 2352 was 2352, checked in by sga, 11 years ago

NEMO branch nemo_v3_3_beta
Add NOCS tools based on SCRIP package for creating weights for interpolation on the fly
These now should build with the maketools script in the TOOLS directory using the same
architecture configuration file as the model (hopefully)

File size: 23.8 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     this program is a short driver that tests the remappings using
4!     a simple analytic field.  the results are written in netCDF
5!     format.
6!
7!-----------------------------------------------------------------------
8!
9!     CVS:$Id: scrip_test_repeat.f,v 1.3 2000/04/19 21:56:26 pwjones Exp $
10!
11!     Copyright (c) 1997, 1998 the Regents of the University of
12!       California.
13!
14!     This software and ancillary information (herein called software)
15!     called SCRIP is made available under the terms described here. 
16!     The software has been approved for release with associated
17!     LA-CC Number 98-45.
18!
19!     Unless otherwise indicated, this software has been authored
20!     by an employee or employees of the University of California,
21!     operator of the Los Alamos National Laboratory under Contract
22!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
23!     Government has rights to use, reproduce, and distribute this
24!     software.  The public may copy and use this software without
25!     charge, provided that this Notice and any statement of authorship
26!     are reproduced on all copies.  Neither the Government nor the
27!     University makes any warranty, express or implied, or assumes
28!     any liability or responsibility for the use of this software.
29!
30!     If software is modified to produce derivative works, such modified
31!     software should be clearly marked, so as not to confuse it with
32!     the version available from Los Alamos National Laboratory.
33!
34!***********************************************************************
35
36      program remap_test_repeat
37
38!-----------------------------------------------------------------------
39
40      use kinds_mod    ! defines common data types
41      use constants    ! defines common constants
42      use iounits      ! I/O unit manager
43      use netcdf_mod   ! netcdf I/O stuff
44      use grids        ! module containing grid info
45      use remap_vars   ! module containing remapping info
46
47      implicit none
48
49!-----------------------------------------------------------------------
50!
51!     interface for remap routine
52!
53!-----------------------------------------------------------------------
54
55      interface
56        subroutine remap(dst_array, map_wts, dst_add, src_add, 
57     &                   src_array, src_grad1, src_grad2, src_grad3)
58
59        use kinds_mod
60        use constants
61
62        implicit none
63
64        integer (kind=int_kind), dimension(:), intent(in) ::
65     &       dst_add, src_add 
66
67        real (kind=dbl_kind), dimension(:,:), intent(in) :: map_wts
68
69        real (kind=dbl_kind), dimension(:), intent(in) :: src_array
70
71        real (kind=dbl_kind), dimension(:), intent(in), optional ::
72     &     src_grad1, src_grad2, src_grad3
73
74        real (kind=dbl_kind), dimension(:), intent(inout) :: dst_array
75
76        end subroutine remap
77      end interface
78
79!-----------------------------------------------------------------------
80!
81!     input namelist variables
82!
83!-----------------------------------------------------------------------
84
85      integer (kind=int_kind) ::
86     &           field_choice, ! choice of field to be interpolated
87     &           num_repeats   ! number of times to repeat remappings
88
89      character (char_len) :: 
90     &           interp_file1, ! filename containing remap data (map1)
91     &           interp_file2, ! filename containing remap data (map2)
92     &           output_file   ! filename containing output test data
93
94      namelist /remap_inputs/ field_choice, num_repeats,
95     &                        interp_file1, interp_file2, output_file
96
97!-----------------------------------------------------------------------
98!
99!     local variables
100!
101!-----------------------------------------------------------------------
102
103      character (char_len) :: 
104     &           map_name1,    ! name for mapping from grid1 to grid2
105     &           map_name2     ! name for mapping from grid2 to grid1
106
107      integer (kind=int_kind) ::    ! netCDF ids for files and arrays
108     &        n, ncstat, nc_outfile_id, 
109     &        nc_srcgrdcntrlat_id, nc_srcgrdcntrlon_id,
110     &        nc_dstgrdcntrlat_id, nc_dstgrdcntrlon_id,
111     &        nc_srcgrdrank_id, nc_dstgrdrank_id,
112     &        nc_srcgrdimask_id, nc_dstgrdimask_id,
113     &        nc_srcgrdarea_id, nc_dstgrdarea_id,
114     &        nc_srcgrdfrac_id, nc_dstgrdfrac_id,
115     &        nc_srcarray_id, nc_srcgradlat_id, nc_srcgradlon_id,
116     &        nc_dstarray1_id, nc_dstarray2_id,
117     &        nc_dsterror1_id, nc_dsterror2_id
118
119      integer (kind=int_kind), dimension(:), allocatable ::
120     &        nc_grid1size_id, nc_grid2size_id
121
122!-----------------------------------------------------------------------
123
124      character (char_len) :: 
125     &          dim_name    ! netCDF dimension name
126
127      integer (kind=int_kind) :: i,j,n,
128     &    iunit  ! unit number for namelist file
129
130      integer (kind=int_kind), dimension(:), allocatable ::
131     &    grid1_imask, grid2_imask
132
133      real (kind=dbl_kind) ::
134     &    length            ! length scale for cosine hill test field
135
136      real (kind=dbl_kind), dimension(:), allocatable ::
137     &    grid1_array, 
138     &    grid1_tmp, 
139     &    grad1_lat, 
140     &    grad1_lon, 
141     &    grid2_array, 
142     &    grid2_err,
143     &    grid2_tmp,
144     &    grad2_lat, 
145     &    grad2_lon 
146
147!-----------------------------------------------------------------------
148!
149!     read namelist for file and mapping info
150!
151!-----------------------------------------------------------------------
152
153      call get_unit(iunit)
154      open(iunit, file='repeat_test_in', status='old', form='formatted')
155      read(iunit, nml=remap_inputs)
156      call release_unit(iunit)
157      write(*,nml=remap_inputs)
158
159!-----------------------------------------------------------------------
160!
161!     read remapping data
162!
163!-----------------------------------------------------------------------
164
165      num_maps = 2
166
167      call read_remap(map_name1, map_name2, interp_file1, interp_file2)
168
169!-----------------------------------------------------------------------
170!
171!     allocate arrays
172!
173!-----------------------------------------------------------------------
174
175      allocate (grid1_array    (grid1_size), 
176     &          grid1_tmp      (grid1_size),
177     &          grad1_lat      (grid1_size), 
178     &          grad1_lon      (grid1_size), 
179     &          grid1_imask    (grid1_size),
180     &          grid2_array    (grid2_size), 
181     &          grid2_err      (grid2_size),
182     &          grid2_tmp      (grid2_size),
183     &          grad2_lat      (grid2_size), 
184     &          grad2_lon      (grid2_size), 
185     &          grid2_imask    (grid2_size))
186
187      where (grid1_mask)
188        grid1_imask = 1
189      elsewhere
190        grid1_imask = 0
191      endwhere
192      where (grid2_mask)
193        grid2_imask = 1
194      elsewhere
195        grid2_imask = 0
196      endwhere
197
198!-----------------------------------------------------------------------
199!
200!     setup a NetCDF file for output
201!
202!-----------------------------------------------------------------------
203
204      !***
205      !*** create netCDF dataset
206      !***
207
208      ncstat = nf_create (output_file, NF_CLOBBER, nc_outfile_id)
209      call netcdf_error_handler(ncstat)
210
211      ncstat = nf_put_att_text (nc_outfile_id, NF_GLOBAL, 'title',
212     &                          len_trim(map_name1), map_name1)
213      call netcdf_error_handler(ncstat)
214
215      !***
216      !*** define grid size dimensions
217      !***
218
219      allocate( nc_grid1size_id(grid1_rank),
220     &          nc_grid2size_id(grid2_rank))
221
222      do n=1,grid1_rank
223        write(dim_name,1000) 'grid1_dim',n
224        ncstat = nf_def_dim (nc_outfile_id, dim_name, 
225     &                       grid1_dims(n), nc_grid1size_id(n))
226        call netcdf_error_handler(ncstat)
227      end do
228
229      do n=1,grid2_rank
230        write(dim_name,1000) 'grid2_dim',n
231        ncstat = nf_def_dim (nc_outfile_id, dim_name, 
232     &                       grid2_dims(n), nc_grid2size_id(n))
233        call netcdf_error_handler(ncstat)
234      end do
235 1000 format(a9,i1)
236
237      !***
238      !*** define grid center latitude array
239      !***
240
241      ncstat = nf_def_var (nc_outfile_id, 'src_grid_center_lat', 
242     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
243     &                     nc_srcgrdcntrlat_id)
244      call netcdf_error_handler(ncstat)
245
246      ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdcntrlat_id, 
247     &                          'units', 7, 'radians')
248      call netcdf_error_handler(ncstat)
249
250      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_center_lat', 
251     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
252     &                     nc_dstgrdcntrlat_id)
253      call netcdf_error_handler(ncstat)
254
255      ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdcntrlat_id, 
256     &                          'units', 7, 'radians')
257      call netcdf_error_handler(ncstat)
258
259      !***
260      !*** define grid center longitude array
261      !***
262
263      ncstat = nf_def_var (nc_outfile_id, 'src_grid_center_lon', 
264     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
265     &                     nc_srcgrdcntrlon_id)
266      call netcdf_error_handler(ncstat)
267
268      ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdcntrlon_id, 
269     &                          'units', 7, 'radians')
270      call netcdf_error_handler(ncstat)
271
272      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_center_lon', 
273     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
274     &                     nc_dstgrdcntrlon_id)
275      call netcdf_error_handler(ncstat)
276
277      ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdcntrlon_id, 
278     &                          'units', 7, 'radians')
279      call netcdf_error_handler(ncstat)
280
281      !***
282      !*** define grid mask
283      !***
284
285      ncstat = nf_def_var (nc_outfile_id, 'src_grid_imask', NF_INT,
286     &               grid1_rank, nc_grid1size_id, nc_srcgrdimask_id)
287      call netcdf_error_handler(ncstat)
288
289      ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdimask_id, 
290     &                          'units', 8, 'unitless')
291      call netcdf_error_handler(ncstat)
292
293      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_imask', NF_INT,
294     &               grid2_rank, nc_grid2size_id, nc_dstgrdimask_id)
295      call netcdf_error_handler(ncstat)
296
297      ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdimask_id, 
298     &                          'units', 8, 'unitless')
299      call netcdf_error_handler(ncstat)
300
301      !***
302      !*** define grid area arrays
303      !***
304
305      ncstat = nf_def_var (nc_outfile_id, 'src_grid_area', 
306     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
307     &                     nc_srcgrdarea_id)
308      call netcdf_error_handler(ncstat)
309
310      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_area', 
311     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
312     &                     nc_dstgrdarea_id)
313      call netcdf_error_handler(ncstat)
314
315      !***
316      !*** define grid fraction arrays
317      !***
318
319      ncstat = nf_def_var (nc_outfile_id, 'src_grid_frac', 
320     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
321     &                     nc_srcgrdfrac_id)
322      call netcdf_error_handler(ncstat)
323
324      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_frac', 
325     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
326     &                     nc_dstgrdfrac_id)
327      call netcdf_error_handler(ncstat)
328
329      !***
330      !*** define source array
331      !***
332
333      ncstat = nf_def_var (nc_outfile_id, 'src_array', 
334     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
335     &                     nc_srcarray_id)
336      call netcdf_error_handler(ncstat)
337
338      !***
339      !*** define gradient arrays
340      !***
341
342      ncstat = nf_def_var (nc_outfile_id, 'src_grad_lat', 
343     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
344     &                     nc_srcgradlat_id)
345      call netcdf_error_handler(ncstat)
346
347      ncstat = nf_def_var (nc_outfile_id, 'src_grad_lon', 
348     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
349     &                     nc_srcgradlon_id)
350      call netcdf_error_handler(ncstat)
351
352      !***
353      !*** define destination arrays
354      !***
355
356      ncstat = nf_def_var (nc_outfile_id, 'dst_array1', 
357     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
358     &                     nc_dstarray1_id)
359      call netcdf_error_handler(ncstat)
360
361      ncstat = nf_def_var (nc_outfile_id, 'dst_array2', 
362     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
363     &                     nc_dstarray2_id)
364      call netcdf_error_handler(ncstat)
365
366      !***
367      !*** define error arrays
368      !***
369
370      ncstat = nf_def_var (nc_outfile_id, 'dst_error1', 
371     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
372     &                     nc_dsterror1_id)
373      call netcdf_error_handler(ncstat)
374
375      ncstat = nf_def_var (nc_outfile_id, 'dst_error2', 
376     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
377     &                     nc_dsterror2_id)
378      call netcdf_error_handler(ncstat)
379
380      !***
381      !*** end definition stage
382      !***
383
384      ncstat = nf_enddef(nc_outfile_id)
385      call netcdf_error_handler(ncstat)
386
387!-----------------------------------------------------------------------
388!
389!     write some grid info
390!
391!-----------------------------------------------------------------------
392
393      !***
394      !*** write grid center latitude array
395      !***
396
397      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdcntrlat_id,
398     &                           grid1_center_lat)
399      call netcdf_error_handler(ncstat)
400
401      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdcntrlat_id,
402     &                           grid2_center_lat)
403      call netcdf_error_handler(ncstat)
404
405      !***
406      !*** write grid center longitude array
407      !***
408
409      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdcntrlon_id,
410     &                           grid1_center_lon)
411      call netcdf_error_handler(ncstat)
412
413      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdcntrlon_id,
414     &                           grid2_center_lon)
415      call netcdf_error_handler(ncstat)
416
417      !***
418      !*** write grid mask
419      !***
420
421      ncstat = nf_put_var_int(nc_outfile_id, nc_srcgrdimask_id,
422     &                        grid1_imask)
423      call netcdf_error_handler(ncstat)
424
425      ncstat = nf_put_var_int(nc_outfile_id, nc_dstgrdimask_id,
426     &                        grid2_imask)
427      call netcdf_error_handler(ncstat)
428
429      !***
430      !*** define grid area arrays
431      !***
432
433      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdarea_id,
434     &                           grid1_area)
435      call netcdf_error_handler(ncstat)
436
437      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdarea_id,
438     &                           grid2_area)
439      call netcdf_error_handler(ncstat)
440
441      !***
442      !*** define grid fraction arrays
443      !***
444
445      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdfrac_id,
446     &                           grid1_frac)
447      call netcdf_error_handler(ncstat)
448
449      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdfrac_id,
450     &                           grid2_frac)
451      call netcdf_error_handler(ncstat)
452
453!-----------------------------------------------------------------------
454!
455!     set up fields for test cases based on user choice
456!
457!-----------------------------------------------------------------------
458
459      select case (field_choice)
460      case(1)  !*** cosine hill at lon=pi and lat=0
461
462        length = 0.1*pi2
463
464        grid1_array = cos(grid1_center_lat)*cos(grid1_center_lon)
465        grid2_array = cos(grid2_center_lat)*cos(grid2_center_lon)
466
467        grid1_tmp = acos(-grid1_array)/length
468        grid2_tmp = acos(-grid2_array)/length
469
470        where (grid1_tmp <= one)
471          grad1_lat   = (pi/length)*sin(pi*grid1_tmp)*
472     &                  sin(grid1_center_lat)*cos(grid1_center_lon)/
473     &                  sqrt(one-grid1_array**2)
474          grad1_lon   = (pi/length)*sin(pi*grid1_tmp)*
475     &                  sin(grid1_center_lon)/
476     &                  sqrt(one-grid1_array**2)
477          grid1_array = two + cos(pi*grid1_tmp)
478        elsewhere
479          grid1_array = one
480          grad1_lat   = zero
481          grad1_lon   = zero
482        endwhere
483       
484        where (grid2_tmp <= one)
485          grad2_lat   = (pi/length)*sin(pi*grid2_tmp)*
486     &                  sin(grid2_center_lat)*cos(grid2_center_lon)/
487     &                  sqrt(one-grid2_array**2)
488          grad2_lon   = (pi/length)*sin(pi*grid2_tmp)*
489     &                  sin(grid2_center_lon)/
490     &                  sqrt(one-grid2_array**2)
491          grid2_array = two + cos(pi*grid2_tmp)
492        elsewhere
493          grid2_array = one
494          grad2_lat   = zero
495          grad2_lon   = zero
496        endwhere
497       
498        where (.not. grid1_mask)
499          grid1_array = zero
500          grad1_lat   = zero
501          grad1_lon   = zero
502        end where
503
504        where (.not. grid2_mask)
505          grid2_array = zero
506          grad2_lat   = zero
507          grad2_lon   = zero
508        end where
509
510      case(2)  !*** pseudo-spherical harmonic l=2,m=2
511
512        where (grid1_mask)
513          grid1_array = two + cos(grid1_center_lat)**2*
514     &                    cos(two*grid1_center_lon)
515          grad1_lat   = -sin(two*grid1_center_lat)*
516     &                   cos(two*grid1_center_lon)
517          grad1_lon   = -two*cos(grid1_center_lat)*
518     &                   sin(two*grid1_center_lon)
519        elsewhere
520          grid1_array = zero
521          grad1_lat   = zero
522          grad1_lon   = zero
523        end where
524
525        where (grid2_mask)
526          grid2_array = two + cos(grid2_center_lat)**2*
527     &                        cos(two*grid2_center_lon)
528          grad2_lat   = -sin(two*grid2_center_lat)*
529     &                   cos(two*grid2_center_lon)
530          grad2_lon   = -two*cos(grid2_center_lat)*
531     &                   sin(two*grid2_center_lon)
532        elsewhere
533          grid2_array = zero
534          grad2_lat   = zero
535          grad2_lon   = zero
536        end where
537
538      case(3)  !*** pseudo-spherical harmonic l=32, m=16
539
540        where (grid1_mask)
541          grid1_array = two + sin(two*grid1_center_lat)**16*
542     &                        cos(16.*grid1_center_lon)
543          grad1_lat   = 32.*sin(two*grid1_center_lat)**15*
544     &                      cos(two*grid1_center_lat)*
545     &                      cos(16.*grid1_center_lon)
546          grad1_lon   = -32.*sin(two*grid1_center_lat)**15*
547     &                       sin(grid1_center_lat)*
548     &                   sin(16.*grid1_center_lon)
549        elsewhere
550          grid1_array = zero
551          grad1_lat   = zero
552          grad1_lon   = zero
553        end where
554
555        where (grid2_mask)
556          grid2_array = two + sin(two*grid2_center_lat)**16*
557     &                        cos(16.*grid2_center_lon)
558          grad2_lat   = 32.*sin(two*grid2_center_lat)**15*
559     &                      cos(two*grid2_center_lat)*
560     &                      cos(16.*grid2_center_lon)
561          grad2_lon   = -32.*sin(two*grid2_center_lat)**15*
562     &                       sin(grid2_center_lat)*
563     &                   sin(16.*grid2_center_lon)
564        elsewhere
565          grid2_array = zero
566          grad2_lat   = zero
567          grad2_lon   = zero
568        end where
569
570      case default
571
572        stop 'Bad choice for field to interpolate'
573
574      end select
575
576!-----------------------------------------------------------------------
577!
578!     test repeated first-order maps between grid1 and grid2
579!
580!-----------------------------------------------------------------------
581
582      call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
583     &           grid1_array)
584      do n=1,num_repeats
585        call remap(grid1_tmp, wts_map2, grid1_add_map2, grid2_add_map2,
586     &             grid2_tmp)
587        call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
588     &             grid1_tmp)
589      end do
590
591      where (grid2_frac > .999)
592        grid2_err = (grid2_tmp - grid2_array)/grid2_array
593      elsewhere 
594        grid2_err = zero
595      end where
596
597      print *,'First order mapping from grid1 to grid2:'
598      print *,'----------------------------------------'
599      print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array)
600      print *,'Grid2 min,max: ',minval(grid2_tmp  ),maxval(grid2_tmp  )
601      print *,' Err2 min,max: ',minval(grid2_err),maxval(grid2_err)
602      print *,' Err2    mean: ',sum(abs(grid2_err))/
603     &                          count(grid2_frac > .001)
604
605      !***
606      !*** Conservation Test
607      !***
608
609      print *,'Conservation:'
610      print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
611      print *,'Grid2 Integral = ',sum(grid2_tmp  *grid2_area*grid2_frac)
612
613!-----------------------------------------------------------------------
614!
615!     write results to NetCDF file
616!
617!-----------------------------------------------------------------------
618
619      ncstat = nf_put_var_double(nc_outfile_id, nc_srcarray_id,
620     &                           grid1_array)
621      call netcdf_error_handler(ncstat)
622
623      ncstat = nf_put_var_double(nc_outfile_id, nc_dstarray1_id,
624     &                           grid2_tmp  )
625      call netcdf_error_handler(ncstat)
626
627      ncstat = nf_put_var_double(nc_outfile_id, nc_dsterror1_id,
628     &                           grid2_err)
629      call netcdf_error_handler(ncstat)
630
631!-----------------------------------------------------------------------
632!
633!     test repeated second-order mapppings between grid1 and grid2
634!
635!-----------------------------------------------------------------------
636
637      if (num_wts > 1) then
638
639      call remap(grid2_tmp  , wts_map1, grid2_add_map1, grid1_add_map1,
640     &           grid1_array, src_grad1=grad1_lat,
641     &                        src_grad2=grad1_lon) 
642
643      do n=1,num_repeats
644        call remap(grid1_tmp, wts_map2, grid1_add_map2, grid2_add_map2,
645     &             grid2_tmp, src_grad1=grad2_lat,
646     &                        src_grad2=grad2_lon) 
647
648        call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
649     &             grid1_tmp, src_grad1=grad1_lat,
650     &                        src_grad2=grad1_lon) 
651
652      end do
653
654      where (grid2_frac > .999)
655        grid2_err = (grid2_tmp - grid2_array)/grid2_array
656      elsewhere 
657        grid2_err = zero
658      end where
659
660      print *,'Second order mapping from grid1 to grid2:'
661      print *,'-----------------------------------------'
662      print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array)
663      print *,'Grid2 min,max: ',minval(grid2_tmp  ),maxval(grid2_tmp  )
664      print *,' Err2 min,max: ',minval(grid2_err),maxval(grid2_err)
665      print *,' Err2    mean: ',sum(abs(grid2_err))/
666     &                          count(grid2_frac > .001)
667
668      !***
669      !*** Conservation Test
670      !***
671
672      print *,'Conservation:'
673      print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
674      print *,'Grid2 Integral = ',sum(grid2_tmp  *grid2_area*grid2_frac)
675
676!-----------------------------------------------------------------------
677!
678!     write results to NetCDF file
679!
680!-----------------------------------------------------------------------
681
682      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgradlon_id,
683     &                           grad1_lon)
684      call netcdf_error_handler(ncstat)
685
686      ncstat = nf_put_var_double(nc_outfile_id, nc_dstarray2_id,
687     &                           grid2_tmp  )
688      call netcdf_error_handler(ncstat)
689
690      ncstat = nf_put_var_double(nc_outfile_id, nc_dsterror2_id,
691     &                           grid2_err)
692      call netcdf_error_handler(ncstat)
693
694      endif
695
696!-----------------------------------------------------------------------
697!
698!     close netCDF file
699!
700!-----------------------------------------------------------------------
701
702      ncstat = nf_close(nc_outfile_id)
703      call netcdf_error_handler(ncstat)
704
705!-----------------------------------------------------------------------
706
707      end program remap_test_repeat
708
709!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.