source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/SCRIP1.4/source/scrip_test.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: 31.4 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!     CVS: $Id: scrip_test.f,v 1.6 2000/04/19 21:45:09 pwjones Exp $
8!
9!-----------------------------------------------------------------------
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
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      use remap_mod    ! module containing remapping routines
47      use remap_read   ! routines for reading remap files
48
49      implicit none
50
51!-----------------------------------------------------------------------
52!
53!     input namelist variables
54!
55!-----------------------------------------------------------------------
56
57      integer (kind=int_kind) ::
58     &  field_choice   ! choice of field to be interpolated
59
60      character (char_len) :: 
61     &  interp_file,   ! filename containing remap data (map1)
62     &  output_file    ! filename for test results
63
64      namelist /remap_inputs/ field_choice, interp_file, output_file
65
66!-----------------------------------------------------------------------
67!
68!     local variables
69!
70!-----------------------------------------------------------------------
71
72      character (char_len) :: 
73     &        map_name      ! name for mapping from grid1 to grid2
74
75      integer (kind=int_kind) ::    ! netCDF ids for files and arrays
76     &        ncstat, nc_outfile_id, 
77     &        nc_srcgrdcntrlat_id, nc_srcgrdcntrlon_id,
78     &        nc_dstgrdcntrlat_id, nc_dstgrdcntrlon_id,
79     &        nc_srcgrdrank_id, nc_dstgrdrank_id,
80     &        nc_srcgrdimask_id, nc_dstgrdimask_id,
81     &        nc_srcgrdarea_id, nc_dstgrdarea_id,
82     &        nc_srcgrdfrac_id, nc_dstgrdfrac_id,
83     &        nc_srcarray_id, nc_srcgradlat_id, nc_srcgradlon_id,
84     &        nc_dstarray1_id, nc_dstarray1a_id, nc_dstarray2_id,
85     &        nc_dsterror1_id, nc_dsterror1a_id, nc_dsterror2_id
86
87      integer (kind=int_kind), dimension(:), allocatable ::
88     &        nc_grid1size_id, nc_grid2size_id
89
90!-----------------------------------------------------------------------
91
92      character (char_len) :: 
93     &          dim_name    ! netCDF dimension name
94
95      integer (kind=int_kind) :: i,j,n,imin,imax,idiff,
96     &    ip1,im1,jp1,jm1,nx,ny, ! for computing bicub gradients
97     &    in,is,ie,iw,ine,inw,ise,isw,
98     &    iunit                  ! unit number for namelist file
99
100      integer (kind=int_kind), dimension(:), allocatable ::
101     &    grid1_imask, grid2_imask, grid2_count
102
103      real (kind=dbl_kind) ::
104     &    delew, delns,     ! variables for computing bicub gradients
105     &    length            ! length scale for cosine hill test field
106
107      real (kind=dbl_kind), dimension(:), allocatable ::
108     &    grid1_array, 
109     &    grid1_tmp, 
110     &    grad1_lat, 
111     &    grad1_lon, 
112     &    grad1_latlon, 
113     &    grad1_lat_zero, 
114     &    grad1_lon_zero, 
115     &    grid2_array, 
116     &    grid2_err,
117     &    grid2_tmp
118
119!-----------------------------------------------------------------------
120!
121!     read namelist for file and mapping info
122!
123!-----------------------------------------------------------------------
124
125      call get_unit(iunit)
126      open(iunit, file='scrip_test_in', status='old', form='formatted')
127      read(iunit, nml=remap_inputs)
128      call release_unit(iunit)
129      write(*,nml=remap_inputs)
130
131!-----------------------------------------------------------------------
132!
133!     read remapping data
134!
135!-----------------------------------------------------------------------
136
137      call read_remap(map_name, interp_file)
138
139!-----------------------------------------------------------------------
140!
141!     allocate arrays
142!
143!-----------------------------------------------------------------------
144
145      allocate (grid1_array    (grid1_size), 
146     &          grid1_tmp      (grid1_size),
147     &          grad1_lat      (grid1_size), 
148     &          grad1_lon      (grid1_size), 
149     &          grad1_lat_zero (grid1_size), 
150     &          grad1_lon_zero (grid1_size), 
151     &          grid1_imask    (grid1_size),
152     &          grid2_array    (grid2_size), 
153     &          grid2_err      (grid2_size),
154     &          grid2_tmp      (grid2_size),
155     &          grid2_imask    (grid2_size),
156     &          grid2_count    (grid2_size))
157
158      where (grid1_mask)
159        grid1_imask = 1
160      elsewhere
161        grid1_imask = 0
162      endwhere
163      where (grid2_mask)
164        grid2_imask = 1
165      elsewhere
166        grid2_imask = 0
167      endwhere
168
169!-----------------------------------------------------------------------
170!
171!     setup a NetCDF file for output
172!
173!-----------------------------------------------------------------------
174
175      !***
176      !*** create netCDF dataset
177      !***
178
179      ncstat = nf_create (output_file, NF_CLOBBER, nc_outfile_id)
180      call netcdf_error_handler(ncstat)
181
182      ncstat = nf_put_att_text (nc_outfile_id, NF_GLOBAL, 'title',
183     &                          len_trim(map_name), map_name)
184      call netcdf_error_handler(ncstat)
185
186      !***
187      !*** define grid size dimensions
188      !***
189
190      allocate( nc_grid1size_id(grid1_rank),
191     &          nc_grid2size_id(grid2_rank))
192
193      do n=1,grid1_rank
194        write(dim_name,1000) 'grid1_dim',n
195        ncstat = nf_def_dim (nc_outfile_id, dim_name, 
196     &                       grid1_dims(n), nc_grid1size_id(n))
197        call netcdf_error_handler(ncstat)
198      end do
199
200      do n=1,grid2_rank
201        write(dim_name,1000) 'grid2_dim',n
202        ncstat = nf_def_dim (nc_outfile_id, dim_name, 
203     &                       grid2_dims(n), nc_grid2size_id(n))
204        call netcdf_error_handler(ncstat)
205      end do
206 1000 format(a9,i1)
207
208      !***
209      !*** define grid center latitude array
210      !***
211
212      ncstat = nf_def_var (nc_outfile_id, 'src_grid_center_lat', 
213     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
214     &                     nc_srcgrdcntrlat_id)
215      call netcdf_error_handler(ncstat)
216
217      ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdcntrlat_id, 
218     &                          'units', 7, 'radians')
219      call netcdf_error_handler(ncstat)
220
221      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_center_lat', 
222     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
223     &                     nc_dstgrdcntrlat_id)
224      call netcdf_error_handler(ncstat)
225
226      ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdcntrlat_id, 
227     &                          'units', 7, 'radians')
228      call netcdf_error_handler(ncstat)
229
230      !***
231      !*** define grid center longitude array
232      !***
233
234      ncstat = nf_def_var (nc_outfile_id, 'src_grid_center_lon', 
235     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
236     &                     nc_srcgrdcntrlon_id)
237      call netcdf_error_handler(ncstat)
238
239      ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdcntrlon_id, 
240     &                          'units', 7, 'radians')
241      call netcdf_error_handler(ncstat)
242
243      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_center_lon', 
244     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
245     &                     nc_dstgrdcntrlon_id)
246      call netcdf_error_handler(ncstat)
247
248      ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdcntrlon_id, 
249     &                          'units', 7, 'radians')
250      call netcdf_error_handler(ncstat)
251
252      !***
253      !*** define grid mask
254      !***
255
256      ncstat = nf_def_var (nc_outfile_id, 'src_grid_imask', NF_INT,
257     &               grid1_rank, nc_grid1size_id, nc_srcgrdimask_id)
258      call netcdf_error_handler(ncstat)
259
260      ncstat = nf_put_att_text (nc_outfile_id, nc_srcgrdimask_id, 
261     &                          'units', 8, 'unitless')
262      call netcdf_error_handler(ncstat)
263
264      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_imask', NF_INT,
265     &               grid2_rank, nc_grid2size_id, nc_dstgrdimask_id)
266      call netcdf_error_handler(ncstat)
267
268      ncstat = nf_put_att_text (nc_outfile_id, nc_dstgrdimask_id, 
269     &                          'units', 8, 'unitless')
270      call netcdf_error_handler(ncstat)
271
272      !***
273      !*** define grid area arrays
274      !***
275
276      ncstat = nf_def_var (nc_outfile_id, 'src_grid_area', 
277     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
278     &                     nc_srcgrdarea_id)
279      call netcdf_error_handler(ncstat)
280
281      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_area', 
282     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
283     &                     nc_dstgrdarea_id)
284      call netcdf_error_handler(ncstat)
285
286      !***
287      !*** define grid fraction arrays
288      !***
289
290      ncstat = nf_def_var (nc_outfile_id, 'src_grid_frac', 
291     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
292     &                     nc_srcgrdfrac_id)
293      call netcdf_error_handler(ncstat)
294
295      ncstat = nf_def_var (nc_outfile_id, 'dst_grid_frac', 
296     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
297     &                     nc_dstgrdfrac_id)
298      call netcdf_error_handler(ncstat)
299
300      !***
301      !*** define source array
302      !***
303
304      ncstat = nf_def_var (nc_outfile_id, 'src_array', 
305     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
306     &                     nc_srcarray_id)
307      call netcdf_error_handler(ncstat)
308
309      !***
310      !*** define gradient arrays
311      !***
312
313      ncstat = nf_def_var (nc_outfile_id, 'src_grad_lat', 
314     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
315     &                     nc_srcgradlat_id)
316      call netcdf_error_handler(ncstat)
317
318      ncstat = nf_def_var (nc_outfile_id, 'src_grad_lon', 
319     &                     NF_DOUBLE, grid1_rank, nc_grid1size_id, 
320     &                     nc_srcgradlon_id)
321      call netcdf_error_handler(ncstat)
322
323      !***
324      !*** define destination arrays
325      !***
326
327      ncstat = nf_def_var (nc_outfile_id, 'dst_array1', 
328     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
329     &                     nc_dstarray1_id)
330      call netcdf_error_handler(ncstat)
331
332      ncstat = nf_def_var (nc_outfile_id, 'dst_array1a', 
333     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
334     &                     nc_dstarray1a_id)
335      call netcdf_error_handler(ncstat)
336
337      ncstat = nf_def_var (nc_outfile_id, 'dst_array2', 
338     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
339     &                     nc_dstarray2_id)
340      call netcdf_error_handler(ncstat)
341
342      !***
343      !*** define error arrays
344      !***
345
346      ncstat = nf_def_var (nc_outfile_id, 'dst_error1', 
347     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
348     &                     nc_dsterror1_id)
349      call netcdf_error_handler(ncstat)
350
351      ncstat = nf_def_var (nc_outfile_id, 'dst_error1a', 
352     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
353     &                     nc_dsterror1a_id)
354      call netcdf_error_handler(ncstat)
355
356      ncstat = nf_def_var (nc_outfile_id, 'dst_error2', 
357     &                     NF_DOUBLE, grid2_rank, nc_grid2size_id, 
358     &                     nc_dsterror2_id)
359      call netcdf_error_handler(ncstat)
360
361      !***
362      !*** end definition stage
363      !***
364
365      ncstat = nf_enddef(nc_outfile_id)
366      call netcdf_error_handler(ncstat)
367
368!-----------------------------------------------------------------------
369!
370!     write some grid info
371!
372!-----------------------------------------------------------------------
373
374      !***
375      !*** write grid center latitude array
376      !***
377
378      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdcntrlat_id,
379     &                           grid1_center_lat)
380      call netcdf_error_handler(ncstat)
381
382      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdcntrlat_id,
383     &                           grid2_center_lat)
384      call netcdf_error_handler(ncstat)
385
386      !***
387      !*** write grid center longitude array
388      !***
389
390      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdcntrlon_id,
391     &                           grid1_center_lon)
392      call netcdf_error_handler(ncstat)
393
394      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdcntrlon_id,
395     &                           grid2_center_lon)
396      call netcdf_error_handler(ncstat)
397
398      !***
399      !*** write grid mask
400      !***
401
402      ncstat = nf_put_var_int(nc_outfile_id, nc_srcgrdimask_id,
403     &                        grid1_imask)
404      call netcdf_error_handler(ncstat)
405
406      ncstat = nf_put_var_int(nc_outfile_id, nc_dstgrdimask_id,
407     &                        grid2_imask)
408      call netcdf_error_handler(ncstat)
409
410      !***
411      !*** define grid area arrays
412      !***
413
414      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdarea_id,
415     &                           grid1_area)
416      call netcdf_error_handler(ncstat)
417
418      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdarea_id,
419     &                           grid2_area)
420      call netcdf_error_handler(ncstat)
421
422      !***
423      !*** define grid fraction arrays
424      !***
425
426      ncstat = nf_put_var_double(nc_outfile_id, nc_srcgrdfrac_id,
427     &                           grid1_frac)
428      call netcdf_error_handler(ncstat)
429
430      ncstat = nf_put_var_double(nc_outfile_id, nc_dstgrdfrac_id,
431     &                           grid2_frac)
432      call netcdf_error_handler(ncstat)
433
434!-----------------------------------------------------------------------
435!
436!     set up fields for test cases based on user choice
437!
438!-----------------------------------------------------------------------
439
440      select case (field_choice)
441      case(1)  !*** cosine hill at lon=pi and lat=0
442
443        length = 0.1*pi2
444
445        grid1_array = cos(grid1_center_lat)*cos(grid1_center_lon)
446        grid2_array = cos(grid2_center_lat)*cos(grid2_center_lon)
447
448        grid1_tmp = acos(-grid1_array)/length
449        grid2_tmp = acos(-grid2_array)/length
450
451        where (grid1_tmp <= one)
452          grad1_lat   = (pi/length)*sin(pi*grid1_tmp)*
453     &                  sin(grid1_center_lat)*cos(grid1_center_lon)/
454     &                  sqrt(one-grid1_array**2)
455          grad1_lon   = (pi/length)*sin(pi*grid1_tmp)*
456     &                  sin(grid1_center_lon)/
457     &                  sqrt(one-grid1_array**2)
458          grid1_array = two + cos(pi*grid1_tmp)
459        elsewhere
460          grid1_array = one
461          grad1_lat   = zero
462          grad1_lon   = zero
463        endwhere
464       
465        where (grid2_tmp <= one)
466          grid2_array = two + cos(pi*grid2_tmp)
467        elsewhere
468          grid2_array = one
469        endwhere
470       
471        where (.not. grid1_mask)
472          grid1_array = zero
473          grad1_lat   = zero
474          grad1_lon   = zero
475        end where
476
477        where (grid2_frac < .001) grid2_array = zero
478
479      case(2)  !*** pseudo-spherical harmonic l=2,m=2
480
481        where (grid1_mask)
482          grid1_array = two + cos(grid1_center_lat)**2*
483     &                    cos(two*grid1_center_lon)
484          grad1_lat   = -sin(two*grid1_center_lat)*
485     &                   cos(two*grid1_center_lon)
486          grad1_lon   = -two*cos(grid1_center_lat)*
487     &                   sin(two*grid1_center_lon)
488        elsewhere
489          grid1_array = zero
490          grad1_lat   = zero
491          grad1_lon   = zero
492        end where
493
494        where (grid2_frac > .001) 
495          grid2_array = two + cos(grid2_center_lat)**2*
496     &                        cos(two*grid2_center_lon)
497        elsewhere
498          grid2_array = zero
499        end where
500
501      case(3)  !*** pseudo-spherical harmonic l=32, m=16
502
503        where (grid1_mask)
504          grid1_array = two + sin(two*grid1_center_lat)**16*
505     &                        cos(16.*grid1_center_lon)
506          grad1_lat   = 32.*sin(two*grid1_center_lat)**15*
507     &                      cos(two*grid1_center_lat)*
508     &                      cos(16.*grid1_center_lon)
509          grad1_lon   = -32.*sin(two*grid1_center_lat)**15*
510     &                       sin(grid1_center_lat)*
511     &                   sin(16.*grid1_center_lon)
512        elsewhere
513          grid1_array = zero
514          grad1_lat   = zero
515          grad1_lon   = zero
516        end where
517
518        where (grid2_frac > .001) 
519          grid2_array = two + sin(two*grid2_center_lat)**16*
520     &                        cos(16.*grid2_center_lon)
521        elsewhere
522          grid2_array = zero
523        end where
524
525      case default
526
527        stop 'Bad choice for field to interpolate'
528
529      end select
530
531!-----------------------------------------------------------------------
532!
533!     if bicubic, we need 3 gradients in logical space
534!
535!-----------------------------------------------------------------------
536
537      if (map_type == map_type_bicubic) then
538        allocate (grad1_latlon (grid1_size)) 
539
540        nx = grid1_dims(1)
541        ny = grid1_dims(2)
542
543        do n=1,grid1_size
544
545          grad1_lat(n) = zero
546          grad1_lon(n) = zero
547          grad1_latlon(n) = zero
548
549          if (grid1_mask(n)) then
550
551            delew = half
552            delns = half
553
554            j = (n-1)/nx + 1
555            i = n - (j-1)*nx
556
557            ip1 = i+1
558            im1 = i-1
559            jp1 = j+1
560            jm1 = j-1
561
562            if (ip1 > nx) ip1 = ip1 - nx
563            if (im1 < 1 ) im1 = nx
564            if (jp1 > ny) then
565              jp1 = j
566              delns = one
567            endif
568            if (jm1 < 1 ) then
569              jm1 = j
570              delns = one
571            endif
572
573            in  = (jp1-1)*nx + i
574            is  = (jm1-1)*nx + i
575            ie  = (j  -1)*nx + ip1
576            iw  = (j  -1)*nx + im1
577
578            ine = (jp1-1)*nx + ip1
579            inw = (jp1-1)*nx + im1
580            ise = (jm1-1)*nx + ip1
581            isw = (jm1-1)*nx + im1
582
583            !*** compute i-gradient
584
585            if (.not. grid1_mask(ie)) then
586              ie = n
587              delew = one
588            endif
589            if (.not. grid1_mask(iw)) then
590              iw = n
591              delew = one
592            endif
593 
594            grad1_lat(n) = delew*(grid1_array(ie) - grid1_array(iw))
595
596            !*** compute j-gradient
597
598            if (.not. grid1_mask(in)) then
599              in = n
600              delns = one
601            endif
602            if (.not. grid1_mask(is)) then
603              is = n
604              delns = one
605            endif
606 
607            grad1_lon(n) = delns*(grid1_array(in) - grid1_array(is))
608
609            !*** compute ij-gradient
610
611            delew = half
612            if (jp1 == j .or. jm1 == j) then
613              delns = one
614            else
615              delns = half
616            endif
617
618            if (.not. grid1_mask(ine)) then
619              if (in /= n) then
620                ine = in
621                delew = one
622              else if (ie /= n) then
623                ine = ie
624                inw = iw
625                if (inw == n) delew = one
626                delns = one
627              else
628                ine = n
629                inw = iw
630                delew = one
631                delns = one
632              endif
633            endif
634
635            if (.not. grid1_mask(inw)) then
636              if (in /= n) then
637                inw = in
638                delew = one
639              else if (iw /= n) then
640                inw = iw
641                ine = ie
642                if (ie == n) delew = one
643                delns = one
644              else
645                inw = n
646                ine = ie
647                delew = one
648                delns = one
649              endif
650            endif
651
652            grad1_lat_zero(n) = delew*(grid1_array(ine) -
653     &                                 grid1_array(inw))
654
655            if (.not. grid1_mask(ise)) then
656              if (is /= n) then
657                ise = is
658                delew = one
659              else if (ie /= n) then
660                ise = ie
661                isw = iw
662                if (isw == n) delew = one
663                delns = one
664              else
665                ise = n
666                isw = iw
667                delew = one
668                delns = one
669              endif
670            endif
671
672            if (.not. grid1_mask(isw)) then
673              if (is /= n) then
674                isw = is
675                delew = one
676              else if (iw /= n) then
677                isw = iw
678                ise = ie
679                if (ie == n) delew = one
680                delns = one
681              else
682                isw = n
683                ise = ie
684                delew = one
685                delns = one
686              endif
687            endif
688
689            grad1_lon_zero(n) = delew*(grid1_array(ise) -
690     &                                 grid1_array(isw))
691
692            grad1_latlon(n) = delns*(grad1_lat_zero(n) -
693     &                               grad1_lon_zero(n))
694
695          endif
696        enddo
697      endif
698
699!-----------------------------------------------------------------------
700!
701!     test a first-order map from grid1 to grid2
702!
703!-----------------------------------------------------------------------
704
705      grad1_lat_zero = zero
706      grad1_lon_zero = zero
707
708      if (map_type /= map_type_bicubic) then
709        call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
710     &             grid1_array)
711      else
712        call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
713     &             grid1_array, src_grad1=grad1_lat,
714     &                          src_grad2=grad1_lon,
715     &                          src_grad3=grad1_latlon)
716      endif
717
718      if (map_type == map_type_conserv) then
719        select case (norm_opt)
720        case (norm_opt_none)
721          grid2_err = grid2_frac*grid2_area
722          where (grid2_err /= zero)
723            grid2_tmp = grid2_tmp/grid2_err
724          else where
725            grid2_tmp = zero
726          end where
727        case (norm_opt_frcarea)
728        case (norm_opt_dstarea)
729          where (grid2_frac /= zero)
730            grid2_tmp = grid2_tmp/grid2_frac
731          else where
732            grid2_tmp = zero
733          end where
734        end select
735      end if
736
737      where (grid2_frac > .999)
738        grid2_err = (grid2_tmp - grid2_array)/grid2_array
739      elsewhere 
740        grid2_err = zero
741      end where
742
743      print *,'First order mapping from grid1 to grid2:'
744      print *,'----------------------------------------'
745      print *,'Grid1 min,max: ',minval(grid1_array),maxval(grid1_array)
746      print *,'Grid2 min,max: ',minval(grid2_tmp  ),maxval(grid2_tmp  )
747      print *,' Err2 min,max: ',minval(grid2_err),maxval(grid2_err)
748      print *,' Err2    mean: ',sum(abs(grid2_err))/
749     &                          count(grid2_frac > .999)
750
751      !***
752      !*** Conservation Test
753      !***
754
755      print *,'Conservation:'
756      print *,'Grid1 Integral = ',sum(grid1_array*grid1_area*grid1_frac)
757      print *,'Grid2 Integral = ',sum(grid2_tmp  *grid2_area*grid2_frac)
758
759!-----------------------------------------------------------------------
760!
761!     write results to NetCDF file
762!
763!-----------------------------------------------------------------------
764
765      ncstat = nf_put_var_double(nc_outfile_id, nc_srcarray_id,
766     &                           grid1_array)
767      call netcdf_error_handler(ncstat)
768
769      ncstat = nf_put_var_double(nc_outfile_id, nc_dstarray1_id,
770     &                           grid2_tmp  )
771      call netcdf_error_handler(ncstat)
772
773      ncstat = nf_put_var_double(nc_outfile_id, nc_dsterror1_id,
774     &                           grid2_err)
775      call netcdf_error_handler(ncstat)
776
777!-----------------------------------------------------------------------
778!
779!     for conservative mappings:
780!     test a second-order map from grid1 to grid2 with only lat grads
781!
782!-----------------------------------------------------------------------
783
784      if (map_type == map_type_conserv) then
785
786        call remap(grid2_tmp, wts_map1, grid2_add_map1, grid1_add_map1,
787     &             grid1_array, src_grad1=grad1_lat,
788     &                          src_grad2=grad1_lon_zero) 
789
790        select case (norm_opt)
791        case (norm_opt_none)
792          grid2_err = grid2_frac*grid2_area
793          where (grid2_err /= zero)
794            grid2_tmp = grid2_tmp/grid2_err
795          else where
796            grid2_tmp = zero
797          end where
798        case (norm_opt_frcarea)
799        case (norm_opt_dstarea)
800          where (grid2_frac /= zero)
801            grid2_tmp = grid2_tmp/grid2_frac
802          else where
803            grid2_tmp = zero
804          end where
805        end select
806
807        where (grid2_frac > .999)
808          grid2_err = (grid2_tmp - grid2_array)/grid2_array
809        elsewhere 
810          grid2_err = zero
811        end where
812
813        print *,'Second order mapping from grid1 to grid2 (lat only):'
814        print *,'----------------------------------------'
815        print *,'Grid1 min,max: ',minval(grid1_array),
816     &                            maxval(grid1_array)
817        print *,'Grid2 min,max: ',minval(grid2_tmp  ),
818     &                            maxval(grid2_tmp  )
819        print *,' Err2 min,max: ',minval(grid2_err),maxval(grid2_err)
820        print *,' Err2    mean: ',sum(abs(grid2_err))/
821     &                            count(grid2_frac > .999)
822
823        !***
824        !*** Conservation Test
825        !***
826
827        print *,'Conservation:'
828        print *,'Grid1 Integral = ',
829     &          sum(grid1_array*grid1_area*grid1_frac)
830        print *,'Grid2 Integral = ',
831     &          sum(grid2_tmp  *grid2_area*grid2_frac)
832
833!-----------------------------------------------------------------------
834!
835!       write results to NetCDF file
836!
837!-----------------------------------------------------------------------
838
839        ncstat = nf_put_var_double(nc_outfile_id, nc_srcgradlat_id,
840     &                             grad1_lat)
841        call netcdf_error_handler(ncstat)
842
843        ncstat = nf_put_var_double(nc_outfile_id, nc_dstarray1a_id,
844     &                             grid2_tmp  )
845        call netcdf_error_handler(ncstat)
846
847        ncstat = nf_put_var_double(nc_outfile_id, nc_dsterror1a_id,
848     &                             grid2_err)
849        call netcdf_error_handler(ncstat)
850
851!-----------------------------------------------------------------------
852!
853!     for conservative mappings:
854!     test a second-order map from grid1 to grid2
855!
856!-----------------------------------------------------------------------
857
858        call remap(grid2_tmp,wts_map1,grid2_add_map1,grid1_add_map1,
859     &             grid1_array, src_grad1=grad1_lat,
860     &                          src_grad2=grad1_lon) 
861
862        select case (norm_opt)
863        case (norm_opt_none)
864          grid2_err = grid2_frac*grid2_area
865          where (grid2_err /= zero)
866            grid2_tmp = grid2_tmp/grid2_err
867          else where
868            grid2_tmp = zero
869          end where
870        case (norm_opt_frcarea)
871        case (norm_opt_dstarea)
872          where (grid2_frac /= zero)
873            grid2_tmp = grid2_tmp/grid2_frac
874          else where
875            grid2_tmp = zero
876          end where
877        end select
878
879        where (grid2_frac > .999)
880          grid2_err = (grid2_tmp - grid2_array)/grid2_array
881        elsewhere 
882          grid2_err = zero
883        end where
884
885        print *,'Second order mapping from grid1 to grid2:'
886        print *,'-----------------------------------------'
887        print *,'Grid1 min,max: ',minval(grid1_array),
888     &                            maxval(grid1_array)
889        print *,'Grid2 min,max: ',minval(grid2_tmp  ),
890     &                            maxval(grid2_tmp  )
891        print *,' Err2 min,max: ',minval(grid2_err),maxval(grid2_err)
892        print *,' Err2    mean: ',sum(abs(grid2_err))/
893     &                            count(grid2_frac > .999)
894
895        !***
896        !*** Conservation Test
897        !***
898
899        print *,'Conservation:'
900        print *,'Grid1 Integral = ',
901     &           sum(grid1_array*grid1_area*grid1_frac)
902        print *,'Grid2 Integral = ',
903     &           sum(grid2_tmp  *grid2_area*grid2_frac)
904
905!-----------------------------------------------------------------------
906!
907!       write results to NetCDF file
908!
909!-----------------------------------------------------------------------
910
911        ncstat = nf_put_var_double(nc_outfile_id, nc_srcgradlon_id,
912     &                             grad1_lon)
913        call netcdf_error_handler(ncstat)
914
915        ncstat = nf_put_var_double(nc_outfile_id, nc_dstarray2_id,
916     &                             grid2_tmp  )
917        call netcdf_error_handler(ncstat)
918
919        ncstat = nf_put_var_double(nc_outfile_id, nc_dsterror2_id,
920     &                             grid2_err)
921        call netcdf_error_handler(ncstat)
922
923      endif
924
925!-----------------------------------------------------------------------
926!
927!     close netCDF file
928!
929!-----------------------------------------------------------------------
930
931      ncstat = nf_close(nc_outfile_id)
932      call netcdf_error_handler(ncstat)
933
934!-----------------------------------------------------------------------
935!
936!     calculate some statistics
937!
938!-----------------------------------------------------------------------
939
940      grid2_count = zero
941      grid2_tmp   = zero
942      grid2_err   = zero
943
944      print *,'number of sparse matrix entries ',num_links_map1
945      do n=1,num_links_map1
946        grid2_count(grid2_add_map1(n)) = 
947     &  grid2_count(grid2_add_map1(n)) + 1
948        if (wts_map1(1,n) > one .or. wts_map1(1,n) < zero) then
949          grid2_tmp(grid2_add_map1(n)) = 
950     &    grid2_tmp(grid2_add_map1(n)) + 1
951          grid2_err(grid2_add_map1(n)) = max(abs(wts_map1(1,n)),
952     &    grid2_err(grid2_add_map1(n)) )
953        endif
954      end do
955
956      do n=1,grid2_size
957        if (grid2_tmp(n) > zero) print *,n,grid2_err(n)
958      end do
959
960      imin = minval(grid2_count, mask=(grid2_count > 0))
961      imax = maxval(grid2_count)
962      idiff =  (imax - imin)/10 + 1
963      print *,'total number of dest cells ',grid2_size
964      print *,'number of cells participating in remap ',
965     &   count(grid2_count > zero)
966      print *,'min no of entries/row = ',imin
967      print *,'max no of entries/row = ',imax
968
969      imax = imin + idiff
970      do n=1,10
971        print *,'num of rows with entries between ',imin,' - ',imax-1,
972     &     count(grid2_count >= imin .and. grid2_count < imax)
973        imin = imin + idiff
974        imax = imax + idiff
975      end do
976
977!-----------------------------------------------------------------------
978
979      end program remap_test
980
981!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.