source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/TOOLS/WEIGHTS/src/remap_write.f90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File size: 56.9 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     This module contains routines for writing the remapping data to
4!     a file.  Before writing the data for each mapping, the links are
5!     sorted by destination grid address.
6!
7!-----------------------------------------------------------------------
8!
9!     CVS:$Id$
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      module remap_write
37
38!-----------------------------------------------------------------------
39
40      use kinds_mod     ! defines common data types
41      use constants     ! defines common scalar constants
42      use grids         ! module containing grid information
43      use remap_vars    ! module containing remap information
44      use netcdf_mod    ! module with netCDF stuff
45
46      implicit none
47
48!-----------------------------------------------------------------------
49!
50!     module variables
51!
52!-----------------------------------------------------------------------
53
54      character(char_len), private ::  &
55         map_method        & ! character string for map_type
56      ,  normalize_opt     & ! character string for normalization option
57      ,  history           & ! character string for history information
58      ,  convention       ! character string for output convention
59
60      character(8), private ::  &
61         cdate            ! character date string
62
63      integer (kind=int_kind), dimension(:), allocatable, private :: &
64         src_mask_int      & ! integer masks to determine
65      ,  dst_mask_int     ! cells that participate in map
66
67!-----------------------------------------------------------------------
68!
69!     various netCDF identifiers used by output routines
70!
71!-----------------------------------------------------------------------
72
73      integer (kind=int_kind), private :: &
74         ncstat                & ! error flag for netCDF calls
75      ,  nc_file_id            & ! id for netCDF file
76      ,  nc_srcgrdsize_id      & ! id for source grid size
77      ,  nc_dstgrdsize_id      & ! id for destination grid size
78      ,  nc_srcgrdcorn_id      & ! id for number of source grid corners
79      ,  nc_dstgrdcorn_id      & ! id for number of dest grid corners
80      ,  nc_srcgrdrank_id      & ! id for source grid rank
81      ,  nc_dstgrdrank_id      & ! id for dest grid rank
82      ,  nc_numlinks_id        & ! id for number of links in mapping
83      ,  nc_numwgts_id         & ! id for number of weights for mapping
84      ,  nc_srcgrddims_id      & ! id for source grid dimensions
85      ,  nc_dstgrddims_id      & ! id for dest grid dimensions
86      ,  nc_srcgrdcntrlat_id   & ! id for source grid center latitude
87      ,  nc_dstgrdcntrlat_id   & ! id for dest grid center latitude
88      ,  nc_srcgrdcntrlon_id   & ! id for source grid center longitude
89      ,  nc_dstgrdcntrlon_id   & ! id for dest grid center longitude
90      ,  nc_srcgrdimask_id     & ! id for source grid mask
91      ,  nc_dstgrdimask_id     & ! id for dest grid mask
92      ,  nc_srcgrdcrnrlat_id   & ! id for latitude of source grid corners
93      ,  nc_srcgrdcrnrlon_id   & ! id for longitude of source grid corners
94      ,  nc_dstgrdcrnrlat_id   & ! id for latitude of dest grid corners
95      ,  nc_dstgrdcrnrlon_id   & ! id for longitude of dest grid corners
96      ,  nc_srcgrdarea_id      & ! id for area of source grid cells
97      ,  nc_dstgrdarea_id      & ! id for area of dest grid cells
98      ,  nc_srcgrdfrac_id      & ! id for area fraction on source grid
99      ,  nc_dstgrdfrac_id      & ! id for area fraction on dest grid
100      ,  nc_srcadd_id          & ! id for map source address
101      ,  nc_dstadd_id          & ! id for map destination address
102      ,  nc_rmpmatrix_id      ! id for remapping matrix
103
104      integer (kind=int_kind), dimension(2), private :: &
105         nc_dims2_id  ! netCDF ids for 2d array dims
106
107!***********************************************************************
108
109      contains
110
111!***********************************************************************
112
113      subroutine write_remap(map1_name, map2_name,  &
114                             interp_file1, interp_file2, output_opt)
115
116!-----------------------------------------------------------------------
117!
118!     calls correct output routine based on output format choice
119!
120!-----------------------------------------------------------------------
121
122!-----------------------------------------------------------------------
123!
124!     input variables
125!
126!-----------------------------------------------------------------------
127
128      character(char_len), intent(in) :: &
129                  map1_name,     & ! name for mapping grid1 to grid2
130                  map2_name,     & ! name for mapping grid2 to grid1
131                  interp_file1,  & ! filename for map1 remap data
132                  interp_file2,  & ! filename for map2 remap data
133                  output_opt    ! option for output conventions
134
135!-----------------------------------------------------------------------
136!
137!     local variables
138!
139!-----------------------------------------------------------------------
140
141!-----------------------------------------------------------------------
142!
143!     define some common variables to be used in all routines
144!
145!-----------------------------------------------------------------------
146
147      select case(norm_opt)
148      case (norm_opt_none)
149        normalize_opt = 'none'
150      case (norm_opt_frcarea)
151        normalize_opt = 'fracarea'
152      case (norm_opt_dstarea)
153        normalize_opt = 'destarea'
154      end select
155
156      select case(map_type)
157      case(map_type_conserv)
158        map_method = 'Conservative remapping'
159      case(map_type_bilinear)
160        map_method = 'Bilinear remapping'
161      case(map_type_distwgt)
162        map_method = 'Distance weighted avg of nearest neighbors'
163      case(map_type_bicubic)
164        map_method = 'Bicubic remapping'
165      case default
166        stop 'Invalid Map Type'
167      end select
168
169      call date_and_time(date=cdate)
170      write (history,1000) cdate(5:6),cdate(7:8),cdate(1:4)
171 1000 format('Created: ',a2,'-',a2,'-',a4)
172
173!-----------------------------------------------------------------------
174!
175!     sort address and weight arrays
176!
177!-----------------------------------------------------------------------
178
179      call sort_add(grid2_add_map1, grid1_add_map1, wts_map1)
180      if (num_maps > 1) then
181        call sort_add(grid1_add_map2, grid2_add_map2, wts_map2)
182      endif
183
184!-----------------------------------------------------------------------
185!
186!     call appropriate output routine
187!
188!-----------------------------------------------------------------------
189
190      select case(output_opt)
191      case ('scrip')
192        call write_remap_scrip(map1_name, interp_file1, 1)
193      case ('ncar-csm')
194        call write_remap_csm  (map1_name, interp_file1, 1)
195      case default
196        stop 'unknown output file convention'
197      end select
198
199!-----------------------------------------------------------------------
200!
201!     call appropriate output routine for second mapping if required
202!
203!-----------------------------------------------------------------------
204
205      if (num_maps > 1) then
206        select case(output_opt)
207        case ('scrip')
208          call write_remap_scrip(map2_name, interp_file2, 2)
209        case ('ncar-csm')
210          call write_remap_csm  (map2_name, interp_file2, 2)
211        case default
212          stop 'unknown output file convention'
213        end select
214      endif
215
216!-----------------------------------------------------------------------
217
218      end subroutine write_remap
219
220!***********************************************************************
221
222      subroutine write_remap_scrip(map_name, interp_file, direction)
223
224!-----------------------------------------------------------------------
225!
226!     writes remap data to a netCDF file using SCRIP conventions
227!
228!-----------------------------------------------------------------------
229
230!-----------------------------------------------------------------------
231!
232!     input variables
233!
234!-----------------------------------------------------------------------
235
236      character(char_len), intent(in) :: &
237                  map_name      & ! name for mapping
238      ,           interp_file  ! filename for remap data
239
240      integer (kind=int_kind), intent(in) :: &
241        direction              ! direction of map (1=grid1 to grid2
242                               !                   2=grid2 to grid1)
243
244!-----------------------------------------------------------------------
245!
246!     local variables
247!
248!-----------------------------------------------------------------------
249
250      character(char_len) :: &
251        grid1_ctmp         & ! character temp for grid1 names
252      , grid2_ctmp        ! character temp for grid2 names
253
254      integer (kind=int_kind) :: &
255        itmp1              & ! integer temp
256      , itmp2              & ! integer temp
257      , itmp3              & ! integer temp
258      , itmp4             ! integer temp
259
260!-----------------------------------------------------------------------
261!
262!     create netCDF file for mapping and define some global attributes
263!
264!-----------------------------------------------------------------------
265
266      ncstat = nf_create (interp_file, NF_CLOBBER, nc_file_id)
267      call netcdf_error_handler(ncstat)
268
269      !***
270      !*** map name
271      !***
272      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, 'title', &
273                                len_trim(map_name), map_name)
274      call netcdf_error_handler(ncstat)
275
276      !***
277      !*** normalization option
278      !***
279      ncstat = nf_put_att_text(nc_file_id, NF_GLOBAL, 'normalization', &
280                               len_trim(normalize_opt), normalize_opt)
281      call netcdf_error_handler(ncstat)
282
283      !***
284      !*** map method
285      !***
286      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, 'map_method', &
287                                len_trim(map_method), map_method)
288      call netcdf_error_handler(ncstat)
289
290      !***
291      !*** history
292      !***
293      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, 'history', &
294                                len_trim(history), history)
295      call netcdf_error_handler(ncstat)
296
297      !***
298      !*** file convention
299      !***
300      convention = 'SCRIP'
301      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, 'conventions', &
302                                len_trim(convention), convention)
303      call netcdf_error_handler(ncstat)
304
305      !***
306      !*** source and destination grid names
307      !***
308
309      if (direction == 1) then
310        grid1_ctmp = 'source_grid'
311        grid2_ctmp = 'dest_grid'
312      else
313        grid1_ctmp = 'dest_grid'
314        grid2_ctmp = 'source_grid'
315      endif
316
317      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, trim(grid1_ctmp), &
318                                len_trim(grid1_name), grid1_name)
319      call netcdf_error_handler(ncstat)
320
321      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, trim(grid2_ctmp), &
322                                len_trim(grid2_name), grid2_name)
323      call netcdf_error_handler(ncstat)
324
325!-----------------------------------------------------------------------
326!
327!     prepare netCDF dimension info
328!
329!-----------------------------------------------------------------------
330
331      !***
332      !*** define grid size dimensions
333      !***
334
335      if (direction == 1) then
336        itmp1 = grid1_size
337        itmp2 = grid2_size
338      else
339        itmp1 = grid2_size
340        itmp2 = grid1_size
341      endif
342
343      ncstat = nf_def_dim (nc_file_id, 'src_grid_size', itmp1,  &
344                           nc_srcgrdsize_id)
345      call netcdf_error_handler(ncstat)
346
347      ncstat = nf_def_dim (nc_file_id, 'dst_grid_size', itmp2,  &
348                           nc_dstgrdsize_id)
349      call netcdf_error_handler(ncstat)
350
351      !***
352      !*** define grid corner dimension
353      !***
354
355      if (direction == 1) then
356        itmp1 = grid1_corners
357        itmp2 = grid2_corners
358      else
359        itmp1 = grid2_corners
360        itmp2 = grid1_corners
361      endif
362
363      ncstat = nf_def_dim (nc_file_id, 'src_grid_corners',  &
364                           itmp1, nc_srcgrdcorn_id)
365      call netcdf_error_handler(ncstat)
366
367      ncstat = nf_def_dim (nc_file_id, 'dst_grid_corners',  &
368                           itmp2, nc_dstgrdcorn_id)
369      call netcdf_error_handler(ncstat)
370
371      !***
372      !*** define grid rank dimension
373      !***
374
375      if (direction == 1) then
376        itmp1 = grid1_rank
377        itmp2 = grid2_rank
378      else
379        itmp1 = grid2_rank
380        itmp2 = grid1_rank
381      endif
382
383      ncstat = nf_def_dim (nc_file_id, 'src_grid_rank',  &
384                           itmp1, nc_srcgrdrank_id)
385      call netcdf_error_handler(ncstat)
386
387      ncstat = nf_def_dim (nc_file_id, 'dst_grid_rank',  &
388                           itmp2, nc_dstgrdrank_id)
389      call netcdf_error_handler(ncstat)
390
391      !***
392      !*** define map size dimensions
393      !***
394
395      if (direction == 1) then
396        itmp1 = num_links_map1
397      else
398        itmp1 = num_links_map2
399      endif
400
401      ncstat = nf_def_dim (nc_file_id, 'num_links',  &
402                           itmp1, nc_numlinks_id)
403      call netcdf_error_handler(ncstat)
404
405      ncstat = nf_def_dim (nc_file_id, 'num_wgts',  &
406                           num_wts, nc_numwgts_id)
407      call netcdf_error_handler(ncstat)
408
409      !***
410      !*** define grid dimensions
411      !***
412
413      ncstat = nf_def_var (nc_file_id, 'src_grid_dims', NF_INT, &
414                           1, nc_srcgrdrank_id, nc_srcgrddims_id)
415      call netcdf_error_handler(ncstat)
416
417      ncstat = nf_def_var (nc_file_id, 'dst_grid_dims', NF_INT, &
418                           1, nc_dstgrdrank_id, nc_dstgrddims_id)
419      call netcdf_error_handler(ncstat)
420
421!-----------------------------------------------------------------------
422!
423!     define all arrays for netCDF descriptors
424!
425!-----------------------------------------------------------------------
426
427      !***
428      !*** define grid center latitude array
429      !***
430
431      ncstat = nf_def_var (nc_file_id, 'src_grid_center_lat',  &
432                           NF_DOUBLE, 1, nc_srcgrdsize_id,  &
433                           nc_srcgrdcntrlat_id)
434      call netcdf_error_handler(ncstat)
435
436      ncstat = nf_def_var (nc_file_id, 'dst_grid_center_lat',  &
437                           NF_DOUBLE, 1, nc_dstgrdsize_id,  &
438                           nc_dstgrdcntrlat_id)
439      call netcdf_error_handler(ncstat)
440
441      !***
442      !*** define grid center longitude array
443      !***
444
445      ncstat = nf_def_var (nc_file_id, 'src_grid_center_lon',  &
446                           NF_DOUBLE, 1, nc_srcgrdsize_id,  &
447                           nc_srcgrdcntrlon_id)
448      call netcdf_error_handler(ncstat)
449
450      ncstat = nf_def_var (nc_file_id, 'dst_grid_center_lon',  &
451                           NF_DOUBLE, 1, nc_dstgrdsize_id,  &
452                           nc_dstgrdcntrlon_id)
453      call netcdf_error_handler(ncstat)
454
455      !***
456      !*** define grid corner lat/lon arrays
457      !***
458
459      nc_dims2_id(1) = nc_srcgrdcorn_id
460      nc_dims2_id(2) = nc_srcgrdsize_id
461
462      ncstat = nf_def_var (nc_file_id, 'src_grid_corner_lat',  &
463                           NF_DOUBLE, 2, nc_dims2_id,  &
464                           nc_srcgrdcrnrlat_id)
465      call netcdf_error_handler(ncstat)
466
467      ncstat = nf_def_var (nc_file_id, 'src_grid_corner_lon',  &
468                           NF_DOUBLE, 2, nc_dims2_id,  &
469                           nc_srcgrdcrnrlon_id)
470      call netcdf_error_handler(ncstat)
471
472      nc_dims2_id(1) = nc_dstgrdcorn_id
473      nc_dims2_id(2) = nc_dstgrdsize_id
474
475      ncstat = nf_def_var (nc_file_id, 'dst_grid_corner_lat',  &
476                           NF_DOUBLE, 2, nc_dims2_id,  &
477                           nc_dstgrdcrnrlat_id)
478      call netcdf_error_handler(ncstat)
479
480      ncstat = nf_def_var (nc_file_id, 'dst_grid_corner_lon',  &
481                           NF_DOUBLE, 2, nc_dims2_id,  &
482                           nc_dstgrdcrnrlon_id)
483      call netcdf_error_handler(ncstat)
484
485      !***
486      !*** define units for all coordinate arrays
487      !***
488
489      if (direction == 1) then
490        grid1_ctmp = grid1_units
491        grid2_ctmp = grid2_units
492      else
493        grid1_ctmp = grid2_units
494        grid2_ctmp = grid1_units
495      endif
496
497      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdcntrlat_id,  &
498                                'units', 7, grid1_ctmp)
499      call netcdf_error_handler(ncstat)
500
501      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdcntrlat_id,  &
502                                'units', 7, grid2_ctmp)
503      call netcdf_error_handler(ncstat)
504
505      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdcntrlon_id,  &
506                                'units', 7, grid1_ctmp)
507      call netcdf_error_handler(ncstat)
508
509      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdcntrlon_id,  &
510                                'units', 7, grid2_ctmp)
511      call netcdf_error_handler(ncstat)
512
513      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdcrnrlat_id,  &
514                                'units', 7, grid1_ctmp)
515      call netcdf_error_handler(ncstat)
516
517      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdcrnrlon_id,  &
518                                'units', 7, grid1_ctmp)
519      call netcdf_error_handler(ncstat)
520
521      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdcrnrlat_id,  &
522                                'units', 7, grid2_ctmp)
523      call netcdf_error_handler(ncstat)
524
525      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdcrnrlon_id,  &
526                                'units', 7, grid2_ctmp)
527      call netcdf_error_handler(ncstat)
528
529      !***
530      !*** define grid mask
531      !***
532
533      ncstat = nf_def_var (nc_file_id, 'src_grid_imask', NF_INT, &
534                           1, nc_srcgrdsize_id, nc_srcgrdimask_id)
535      call netcdf_error_handler(ncstat)
536
537      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdimask_id,  &
538                                'units', 8, 'unitless')
539      call netcdf_error_handler(ncstat)
540
541      ncstat = nf_def_var (nc_file_id, 'dst_grid_imask', NF_INT, &
542                           1, nc_dstgrdsize_id, nc_dstgrdimask_id)
543      call netcdf_error_handler(ncstat)
544
545      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdimask_id,  &
546                                'units', 8, 'unitless')
547      call netcdf_error_handler(ncstat)
548
549      !***
550      !*** define grid area arrays
551      !***
552
553      ncstat = nf_def_var (nc_file_id, 'src_grid_area',  &
554                           NF_DOUBLE, 1, nc_srcgrdsize_id,  &
555                           nc_srcgrdarea_id)
556      call netcdf_error_handler(ncstat)
557
558      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdarea_id,  &
559                                'units', 14, 'square radians')
560      call netcdf_error_handler(ncstat)
561
562      ncstat = nf_def_var (nc_file_id, 'dst_grid_area',  &
563                           NF_DOUBLE, 1, nc_dstgrdsize_id,  &
564                           nc_dstgrdarea_id)
565      call netcdf_error_handler(ncstat)
566
567      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdarea_id,  &
568                                'units', 14, 'square radians')
569      call netcdf_error_handler(ncstat)
570
571      !***
572      !*** define grid fraction arrays
573      !***
574
575      ncstat = nf_def_var (nc_file_id, 'src_grid_frac',  &
576                           NF_DOUBLE, 1, nc_srcgrdsize_id,  &
577                           nc_srcgrdfrac_id)
578      call netcdf_error_handler(ncstat)
579
580      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdfrac_id,  &
581                                'units', 8, 'unitless')
582      call netcdf_error_handler(ncstat)
583
584      ncstat = nf_def_var (nc_file_id, 'dst_grid_frac',  &
585                           NF_DOUBLE, 1, nc_dstgrdsize_id,  &
586                           nc_dstgrdfrac_id)
587      call netcdf_error_handler(ncstat)
588
589      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdfrac_id,  &
590                                'units', 8, 'unitless')
591      call netcdf_error_handler(ncstat)
592
593      !***
594      !*** define mapping arrays
595      !***
596
597      ncstat = nf_def_var (nc_file_id, 'src_address',  &
598                           NF_INT, 1, nc_numlinks_id,  &
599                           nc_srcadd_id)
600      call netcdf_error_handler(ncstat)
601
602      ncstat = nf_def_var (nc_file_id, 'dst_address',  &
603                           NF_INT, 1, nc_numlinks_id,  &
604                           nc_dstadd_id)
605      call netcdf_error_handler(ncstat)
606
607      nc_dims2_id(1) = nc_numwgts_id
608      nc_dims2_id(2) = nc_numlinks_id
609
610      ncstat = nf_def_var (nc_file_id, 'remap_matrix',  &
611                           NF_DOUBLE, 2, nc_dims2_id,  &
612                           nc_rmpmatrix_id)
613      call netcdf_error_handler(ncstat)
614
615      !***
616      !*** end definition stage
617      !***
618
619      ncstat = nf_enddef(nc_file_id)
620      call netcdf_error_handler(ncstat)
621
622!-----------------------------------------------------------------------
623!
624!     compute integer masks
625!
626!-----------------------------------------------------------------------
627
628      if (direction == 1) then
629        allocate (src_mask_int(grid1_size), &
630                  dst_mask_int(grid2_size))
631
632        where (grid2_mask)
633          dst_mask_int = 1
634        elsewhere
635          dst_mask_int = 0
636        endwhere
637
638        where (grid1_mask)
639          src_mask_int = 1
640        elsewhere
641          src_mask_int = 0
642        endwhere
643      else
644        allocate (src_mask_int(grid2_size), &
645                  dst_mask_int(grid1_size))
646
647        where (grid1_mask)
648          dst_mask_int = 1
649        elsewhere
650          dst_mask_int = 0
651        endwhere
652
653        where (grid2_mask)
654          src_mask_int = 1
655        elsewhere
656          src_mask_int = 0
657        endwhere
658      endif
659
660!-----------------------------------------------------------------------
661!
662!     change units of lat/lon coordinates if input units different
663!     from radians
664!
665!-----------------------------------------------------------------------
666
667      if (grid1_units(1:7) == 'degrees' .and. direction == 1) then
668        grid1_center_lat = grid1_center_lat/deg2rad
669        grid1_center_lon = grid1_center_lon/deg2rad
670        grid1_corner_lat = grid1_corner_lat/deg2rad
671        grid1_corner_lon = grid1_corner_lon/deg2rad
672      endif
673
674      if (grid2_units(1:7) == 'degrees' .and. direction == 1) then
675        grid2_center_lat = grid2_center_lat/deg2rad
676        grid2_center_lon = grid2_center_lon/deg2rad
677        grid2_corner_lat = grid2_corner_lat/deg2rad
678        grid2_corner_lon = grid2_corner_lon/deg2rad
679      endif
680
681!-----------------------------------------------------------------------
682!
683!     write mapping data
684!
685!-----------------------------------------------------------------------
686
687      if (direction == 1) then
688        itmp1 = nc_srcgrddims_id
689        itmp2 = nc_dstgrddims_id
690      else
691        itmp2 = nc_srcgrddims_id
692        itmp1 = nc_dstgrddims_id
693      endif
694
695      ncstat = nf_put_var_int(nc_file_id, itmp1, grid1_dims)
696      call netcdf_error_handler(ncstat)
697
698      ncstat = nf_put_var_int(nc_file_id, itmp2, grid2_dims)
699      call netcdf_error_handler(ncstat)
700
701      ncstat = nf_put_var_int(nc_file_id, nc_srcgrdimask_id,  &
702                              src_mask_int)
703      call netcdf_error_handler(ncstat)
704
705      ncstat = nf_put_var_int(nc_file_id, nc_dstgrdimask_id, &
706                              dst_mask_int)
707      call netcdf_error_handler(ncstat)
708
709      deallocate(src_mask_int, dst_mask_int)
710
711      if (direction == 1) then
712        itmp1 = nc_srcgrdcntrlat_id
713        itmp2 = nc_srcgrdcntrlon_id
714        itmp3 = nc_srcgrdcrnrlat_id
715        itmp4 = nc_srcgrdcrnrlon_id
716      else
717        itmp1 = nc_dstgrdcntrlat_id
718        itmp2 = nc_dstgrdcntrlon_id
719        itmp3 = nc_dstgrdcrnrlat_id
720        itmp4 = nc_dstgrdcrnrlon_id
721      endif
722
723      ncstat = nf_put_var_double(nc_file_id, itmp1, grid1_center_lat)
724      call netcdf_error_handler(ncstat)
725
726      ncstat = nf_put_var_double(nc_file_id, itmp2, grid1_center_lon)
727      call netcdf_error_handler(ncstat)
728
729      ncstat = nf_put_var_double(nc_file_id, itmp3, grid1_corner_lat)
730      call netcdf_error_handler(ncstat)
731
732      ncstat = nf_put_var_double(nc_file_id, itmp4, grid1_corner_lon)
733      call netcdf_error_handler(ncstat)
734
735      if (direction == 1) then
736        itmp1 = nc_dstgrdcntrlat_id
737        itmp2 = nc_dstgrdcntrlon_id
738        itmp3 = nc_dstgrdcrnrlat_id
739        itmp4 = nc_dstgrdcrnrlon_id
740      else
741        itmp1 = nc_srcgrdcntrlat_id
742        itmp2 = nc_srcgrdcntrlon_id
743        itmp3 = nc_srcgrdcrnrlat_id
744        itmp4 = nc_srcgrdcrnrlon_id
745      endif
746
747      ncstat = nf_put_var_double(nc_file_id, itmp1, grid2_center_lat)
748      call netcdf_error_handler(ncstat)
749
750      ncstat = nf_put_var_double(nc_file_id, itmp2, grid2_center_lon)
751      call netcdf_error_handler(ncstat)
752
753      ncstat = nf_put_var_double(nc_file_id, itmp3, grid2_corner_lat)
754      call netcdf_error_handler(ncstat)
755
756      ncstat = nf_put_var_double(nc_file_id, itmp4, grid2_corner_lon)
757      call netcdf_error_handler(ncstat)
758
759      if (direction == 1) then
760        itmp1 = nc_srcgrdarea_id
761        itmp2 = nc_srcgrdfrac_id
762        itmp3 = nc_dstgrdarea_id
763        itmp4 = nc_dstgrdfrac_id
764      else
765        itmp1 = nc_dstgrdarea_id
766        itmp2 = nc_dstgrdfrac_id
767        itmp3 = nc_srcgrdarea_id
768        itmp4 = nc_srcgrdfrac_id
769      endif
770
771      if (luse_grid1_area) then
772        ncstat = nf_put_var_double(nc_file_id, itmp1, grid1_area_in)
773      else
774        ncstat = nf_put_var_double(nc_file_id, itmp1, grid1_area)
775      endif
776      call netcdf_error_handler(ncstat)
777
778      ncstat = nf_put_var_double(nc_file_id, itmp2, grid1_frac)
779      call netcdf_error_handler(ncstat)
780
781      if (luse_grid2_area) then
782        ncstat = nf_put_var_double(nc_file_id, itmp3, grid2_area_in)
783      else
784        ncstat = nf_put_var_double(nc_file_id, itmp3, grid2_area)
785      endif
786      call netcdf_error_handler(ncstat)
787
788      ncstat = nf_put_var_double(nc_file_id, itmp4, grid2_frac)
789      call netcdf_error_handler(ncstat)
790
791      if (direction == 1) then
792        ncstat = nf_put_var_int(nc_file_id, nc_srcadd_id,  &
793                                grid1_add_map1)
794        call netcdf_error_handler(ncstat)
795
796        ncstat = nf_put_var_int(nc_file_id, nc_dstadd_id,  &
797                                grid2_add_map1)
798        call netcdf_error_handler(ncstat)
799
800        ncstat = nf_put_var_double(nc_file_id, nc_rmpmatrix_id,  &
801                                   wts_map1)
802        call netcdf_error_handler(ncstat)
803      else
804        ncstat = nf_put_var_int(nc_file_id, nc_srcadd_id,  &
805                                grid2_add_map2)
806        call netcdf_error_handler(ncstat)
807
808        ncstat = nf_put_var_int(nc_file_id, nc_dstadd_id,  &
809                                grid1_add_map2)
810        call netcdf_error_handler(ncstat)
811
812        ncstat = nf_put_var_double(nc_file_id, nc_rmpmatrix_id,  &
813                                   wts_map2)
814        call netcdf_error_handler(ncstat)
815      endif
816
817      ncstat = nf_close(nc_file_id)
818      call netcdf_error_handler(ncstat)
819
820!-----------------------------------------------------------------------
821
822      end subroutine write_remap_scrip
823
824!***********************************************************************
825
826      subroutine write_remap_csm(map_name, interp_file, direction)
827
828!-----------------------------------------------------------------------
829!
830!     writes remap data to a netCDF file using NCAR-CSM conventions
831!
832!-----------------------------------------------------------------------
833
834!-----------------------------------------------------------------------
835!
836!     input variables
837!
838!-----------------------------------------------------------------------
839
840      character(char_len), intent(in) :: &
841                  map_name      & ! name for mapping
842      ,           interp_file  ! filename for remap data
843
844      integer (kind=int_kind), intent(in) :: &
845        direction              ! direction of map (1=grid1 to grid2
846                               !                   2=grid2 to grid1)
847
848!-----------------------------------------------------------------------
849!
850!     local variables
851!
852!-----------------------------------------------------------------------
853
854      character(char_len) :: &
855        grid1_ctmp         & ! character temp for grid1 names
856      , grid2_ctmp        ! character temp for grid2 names
857
858      integer (kind=int_kind) :: &
859        itmp1              & ! integer temp
860      , itmp2              & ! integer temp
861      , itmp3              & ! integer temp
862      , itmp4              & ! integer temp
863      , nc_numwgts1_id     & ! extra netCDF id for additional weights
864      , nc_src_isize_id    & ! extra netCDF id for ni_a
865      , nc_src_jsize_id    & ! extra netCDF id for nj_a
866      , nc_dst_isize_id    & ! extra netCDF id for ni_b
867      , nc_dst_jsize_id    & ! extra netCDF id for nj_b
868      , nc_rmpmatrix2_id  ! extra netCDF id for high-order remap matrix
869
870      real (kind=dbl_kind), dimension(:),allocatable :: &
871        wts1              ! CSM wants single array for 1st-order wts
872
873      real (kind=dbl_kind), dimension(:,:),allocatable :: &
874        wts2              ! write remaining weights in different array
875
876!-----------------------------------------------------------------------
877!
878!     create netCDF file for mapping and define some global attributes
879!
880!-----------------------------------------------------------------------
881
882      ncstat = nf_create (interp_file, NF_CLOBBER, nc_file_id)
883      call netcdf_error_handler(ncstat)
884
885      !***
886      !*** map name
887      !***
888      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, 'title', &
889                                len_trim(map_name), map_name)
890      call netcdf_error_handler(ncstat)
891
892      !***
893      !*** normalization option
894      !***
895      ncstat = nf_put_att_text(nc_file_id, NF_GLOBAL, 'normalization', &
896                               len_trim(normalize_opt), normalize_opt)
897      call netcdf_error_handler(ncstat)
898
899      !***
900      !*** map method
901      !***
902      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, 'map_method', &
903                                len_trim(map_method), map_method)
904      call netcdf_error_handler(ncstat)
905
906      !***
907      !*** history
908      !***
909      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, 'history', &
910                                len_trim(history), history)
911      call netcdf_error_handler(ncstat)
912
913      !***
914      !*** file convention
915      !***
916      convention = 'NCAR-CSM'
917      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, 'conventions', &
918                                len_trim(convention), convention)
919      call netcdf_error_handler(ncstat)
920
921      !***
922      !*** source and destination grid names
923      !***
924
925      if (direction == 1) then
926        grid1_ctmp = 'domain_a'
927        grid2_ctmp = 'domain_b'
928      else
929        grid1_ctmp = 'domain_b'
930        grid2_ctmp = 'domain_a'
931      endif
932
933      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, trim(grid1_ctmp), &
934                                len_trim(grid1_name), grid1_name)
935      call netcdf_error_handler(ncstat)
936
937      ncstat = nf_put_att_text (nc_file_id, NF_GLOBAL, trim(grid2_ctmp), &
938                                len_trim(grid2_name), grid2_name)
939      call netcdf_error_handler(ncstat)
940
941!-----------------------------------------------------------------------
942!
943!     prepare netCDF dimension info
944!
945!-----------------------------------------------------------------------
946
947      !***
948      !*** define grid size dimensions
949      !***
950
951      if (direction == 1) then
952        itmp1 = grid1_size
953        itmp2 = grid2_size
954      else
955        itmp1 = grid2_size
956        itmp2 = grid1_size
957      endif
958
959      ncstat = nf_def_dim (nc_file_id, 'n_a', itmp1, nc_srcgrdsize_id)
960      call netcdf_error_handler(ncstat)
961
962      ncstat = nf_def_dim (nc_file_id, 'n_b', itmp2, nc_dstgrdsize_id)
963      call netcdf_error_handler(ncstat)
964
965      !***
966      !*** define grid corner dimension
967      !***
968
969      if (direction == 1) then
970        itmp1 = grid1_corners
971        itmp2 = grid2_corners
972      else
973        itmp1 = grid2_corners
974        itmp2 = grid1_corners
975      endif
976
977      ncstat = nf_def_dim (nc_file_id, 'nv_a', itmp1, nc_srcgrdcorn_id)
978      call netcdf_error_handler(ncstat)
979
980      ncstat = nf_def_dim (nc_file_id, 'nv_b', itmp2, nc_dstgrdcorn_id)
981      call netcdf_error_handler(ncstat)
982
983      !***
984      !*** define grid rank dimension
985      !***
986
987      if (direction == 1) then
988        itmp1 = grid1_rank
989        itmp2 = grid2_rank
990      else
991        itmp1 = grid2_rank
992        itmp2 = grid1_rank
993      endif
994
995      ncstat = nf_def_dim (nc_file_id, 'src_grid_rank',  &
996                           itmp1, nc_srcgrdrank_id)
997      call netcdf_error_handler(ncstat)
998
999      ncstat = nf_def_dim (nc_file_id, 'dst_grid_rank',  &
1000                           itmp2, nc_dstgrdrank_id)
1001      call netcdf_error_handler(ncstat)
1002
1003      !***
1004      !*** define first two dims as if 2-d cartesian domain
1005      !***
1006
1007      if (direction == 1) then
1008        itmp1 = grid1_dims(1)
1009        if (grid1_rank > 1) then
1010          itmp2 = grid1_dims(2)
1011        else
1012          itmp2 = 0
1013        endif
1014        itmp3 = grid2_dims(1)
1015        if (grid2_rank > 1) then
1016          itmp4 = grid2_dims(2)
1017        else
1018          itmp4 = 0
1019        endif
1020      else
1021        itmp1 = grid2_dims(1)
1022        if (grid2_rank > 1) then
1023          itmp2 = grid2_dims(2)
1024        else
1025          itmp2 = 0
1026        endif
1027        itmp3 = grid1_dims(1)
1028        if (grid1_rank > 1) then
1029          itmp4 = grid1_dims(2)
1030        else
1031          itmp4 = 0
1032        endif
1033      endif
1034
1035      ncstat = nf_def_dim (nc_file_id, 'ni_a', itmp1, nc_src_isize_id)
1036      call netcdf_error_handler(ncstat)
1037
1038      ncstat = nf_def_dim (nc_file_id, 'nj_a', itmp2, nc_src_jsize_id)
1039      call netcdf_error_handler(ncstat)
1040
1041      ncstat = nf_def_dim (nc_file_id, 'ni_b', itmp3, nc_dst_isize_id)
1042      call netcdf_error_handler(ncstat)
1043
1044      ncstat = nf_def_dim (nc_file_id, 'nj_b', itmp4, nc_dst_jsize_id)
1045      call netcdf_error_handler(ncstat)
1046
1047      !***
1048      !*** define map size dimensions
1049      !***
1050
1051      if (direction == 1) then
1052        itmp1 = num_links_map1
1053      else
1054        itmp1 = num_links_map2
1055      endif
1056
1057      ncstat = nf_def_dim (nc_file_id, 'n_s', itmp1, nc_numlinks_id)
1058      call netcdf_error_handler(ncstat)
1059
1060      ncstat = nf_def_dim (nc_file_id, 'num_wgts',  &
1061                           num_wts, nc_numwgts_id)
1062      call netcdf_error_handler(ncstat)
1063
1064      if (num_wts > 1) then
1065        ncstat = nf_def_dim (nc_file_id, 'num_wgts1',  &
1066                             num_wts-1, nc_numwgts1_id)
1067        call netcdf_error_handler(ncstat)
1068      endif
1069
1070      !***
1071      !*** define grid dimensions
1072      !***
1073
1074      ncstat = nf_def_var (nc_file_id, 'src_grid_dims', NF_INT, &
1075                           1, nc_srcgrdrank_id, nc_srcgrddims_id)
1076      call netcdf_error_handler(ncstat)
1077
1078      ncstat = nf_def_var (nc_file_id, 'dst_grid_dims', NF_INT, &
1079                           1, nc_dstgrdrank_id, nc_dstgrddims_id)
1080      call netcdf_error_handler(ncstat)
1081
1082!-----------------------------------------------------------------------
1083!
1084!     define all arrays for netCDF descriptors
1085!
1086!-----------------------------------------------------------------------
1087
1088      !***
1089      !*** define grid center latitude array
1090      !***
1091
1092      ncstat = nf_def_var (nc_file_id, 'yc_a', &
1093                           NF_DOUBLE, 1, nc_srcgrdsize_id,  &
1094                           nc_srcgrdcntrlat_id)
1095      call netcdf_error_handler(ncstat)
1096
1097      ncstat = nf_def_var (nc_file_id, 'yc_b',  &
1098                           NF_DOUBLE, 1, nc_dstgrdsize_id,  &
1099                           nc_dstgrdcntrlat_id)
1100      call netcdf_error_handler(ncstat)
1101
1102      !***
1103      !*** define grid center longitude array
1104      !***
1105
1106      ncstat = nf_def_var (nc_file_id, 'xc_a',  &
1107                           NF_DOUBLE, 1, nc_srcgrdsize_id,  &
1108                           nc_srcgrdcntrlon_id)
1109      call netcdf_error_handler(ncstat)
1110
1111      ncstat = nf_def_var (nc_file_id, 'xc_b',  &
1112                           NF_DOUBLE, 1, nc_dstgrdsize_id,  &
1113                           nc_dstgrdcntrlon_id)
1114      call netcdf_error_handler(ncstat)
1115
1116      !***
1117      !*** define grid corner lat/lon arrays
1118      !***
1119
1120      nc_dims2_id(1) = nc_srcgrdcorn_id
1121      nc_dims2_id(2) = nc_srcgrdsize_id
1122
1123      ncstat = nf_def_var (nc_file_id, 'yv_a',  &
1124                           NF_DOUBLE, 2, nc_dims2_id,  &
1125                           nc_srcgrdcrnrlat_id)
1126      call netcdf_error_handler(ncstat)
1127
1128      ncstat = nf_def_var (nc_file_id, 'xv_a',  &
1129                           NF_DOUBLE, 2, nc_dims2_id,  &
1130                           nc_srcgrdcrnrlon_id)
1131      call netcdf_error_handler(ncstat)
1132
1133      nc_dims2_id(1) = nc_dstgrdcorn_id
1134      nc_dims2_id(2) = nc_dstgrdsize_id
1135
1136      ncstat = nf_def_var (nc_file_id, 'yv_b',  &
1137                           NF_DOUBLE, 2, nc_dims2_id,  &
1138                           nc_dstgrdcrnrlat_id)
1139      call netcdf_error_handler(ncstat)
1140
1141      ncstat = nf_def_var (nc_file_id, 'xv_b',  &
1142                           NF_DOUBLE, 2, nc_dims2_id,  &
1143                           nc_dstgrdcrnrlon_id)
1144      call netcdf_error_handler(ncstat)
1145
1146      !***
1147      !*** CSM wants all in degrees
1148      !***
1149
1150      grid1_units = 'degrees'
1151      grid2_units = 'degrees'
1152
1153      if (direction == 1) then
1154        grid1_ctmp = grid1_units
1155        grid2_ctmp = grid2_units
1156      else
1157        grid1_ctmp = grid2_units
1158        grid2_ctmp = grid1_units
1159      endif
1160
1161      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdcntrlat_id,  &
1162                                'units', 7, grid1_ctmp)
1163      call netcdf_error_handler(ncstat)
1164
1165      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdcntrlat_id,  &
1166                                'units', 7, grid2_ctmp)
1167      call netcdf_error_handler(ncstat)
1168
1169      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdcntrlon_id,  &
1170                                'units', 7, grid1_ctmp)
1171      call netcdf_error_handler(ncstat)
1172
1173      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdcntrlon_id,  &
1174                                'units', 7, grid2_ctmp)
1175      call netcdf_error_handler(ncstat)
1176
1177      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdcrnrlat_id,  &
1178                                'units', 7, grid1_ctmp)
1179      call netcdf_error_handler(ncstat)
1180
1181      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdcrnrlon_id,  &
1182                                'units', 7, grid1_ctmp)
1183      call netcdf_error_handler(ncstat)
1184
1185      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdcrnrlat_id,  &
1186                                'units', 7, grid2_ctmp)
1187      call netcdf_error_handler(ncstat)
1188
1189      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdcrnrlon_id,  &
1190                                'units', 7, grid2_ctmp)
1191      call netcdf_error_handler(ncstat)
1192
1193      !***
1194      !*** define grid mask
1195      !***
1196
1197      ncstat = nf_def_var (nc_file_id, 'mask_a', NF_INT, &
1198                           1, nc_srcgrdsize_id, nc_srcgrdimask_id)
1199      call netcdf_error_handler(ncstat)
1200
1201      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdimask_id,  &
1202                                'units', 8, 'unitless')
1203      call netcdf_error_handler(ncstat)
1204
1205      ncstat = nf_def_var (nc_file_id, 'mask_b', NF_INT, &
1206                           1, nc_dstgrdsize_id, nc_dstgrdimask_id)
1207      call netcdf_error_handler(ncstat)
1208
1209      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdimask_id,  &
1210                                'units', 8, 'unitless')
1211      call netcdf_error_handler(ncstat)
1212
1213      !***
1214      !*** define grid area arrays
1215      !***
1216
1217      ncstat = nf_def_var (nc_file_id, 'area_a',  &
1218                           NF_DOUBLE, 1, nc_srcgrdsize_id,  &
1219                           nc_srcgrdarea_id)
1220      call netcdf_error_handler(ncstat)
1221
1222      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdarea_id,  &
1223                                'units', 14, 'square radians')
1224      call netcdf_error_handler(ncstat)
1225
1226      ncstat = nf_def_var (nc_file_id, 'area_b',  &
1227                           NF_DOUBLE, 1, nc_dstgrdsize_id,  &
1228                           nc_dstgrdarea_id)
1229      call netcdf_error_handler(ncstat)
1230
1231      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdarea_id,  &
1232                                'units', 14, 'square radians')
1233      call netcdf_error_handler(ncstat)
1234
1235      !***
1236      !*** define grid fraction arrays
1237      !***
1238
1239      ncstat = nf_def_var (nc_file_id, 'frac_a',  &
1240                           NF_DOUBLE, 1, nc_srcgrdsize_id,  &
1241                           nc_srcgrdfrac_id)
1242      call netcdf_error_handler(ncstat)
1243
1244      ncstat = nf_put_att_text (nc_file_id, nc_srcgrdfrac_id,  &
1245                                'units', 8, 'unitless')
1246      call netcdf_error_handler(ncstat)
1247
1248      ncstat = nf_def_var (nc_file_id, 'frac_b',  &
1249                           NF_DOUBLE, 1, nc_dstgrdsize_id,  &
1250                           nc_dstgrdfrac_id)
1251      call netcdf_error_handler(ncstat)
1252
1253      ncstat = nf_put_att_text (nc_file_id, nc_dstgrdfrac_id,  &
1254                                'units', 8, 'unitless')
1255      call netcdf_error_handler(ncstat)
1256
1257      !***
1258      !*** define mapping arrays
1259      !***
1260
1261      ncstat = nf_def_var (nc_file_id, 'col',  &
1262                           NF_INT, 1, nc_numlinks_id,  &
1263                           nc_srcadd_id)
1264      call netcdf_error_handler(ncstat)
1265
1266      ncstat = nf_def_var (nc_file_id, 'row',  &
1267                           NF_INT, 1, nc_numlinks_id,  &
1268                           nc_dstadd_id)
1269      call netcdf_error_handler(ncstat)
1270
1271      ncstat = nf_def_var (nc_file_id, 'S',  &
1272                           NF_DOUBLE, 1, nc_numlinks_id,  &
1273                           nc_rmpmatrix_id)
1274      call netcdf_error_handler(ncstat)
1275
1276      if (num_wts > 1) then
1277        nc_dims2_id(1) = nc_numwgts1_id
1278        nc_dims2_id(2) = nc_numlinks_id
1279
1280        ncstat = nf_def_var (nc_file_id, 'S2',  &
1281                           NF_DOUBLE, 2, nc_dims2_id,  &
1282                           nc_rmpmatrix2_id)
1283        call netcdf_error_handler(ncstat)
1284      endif
1285
1286      !***
1287      !*** end definition stage
1288      !***
1289
1290      ncstat = nf_enddef(nc_file_id)
1291      call netcdf_error_handler(ncstat)
1292
1293!-----------------------------------------------------------------------
1294!
1295!     compute integer masks
1296!
1297!-----------------------------------------------------------------------
1298
1299      if (direction == 1) then
1300        allocate (src_mask_int(grid1_size), &
1301                  dst_mask_int(grid2_size))
1302
1303        where (grid2_mask)
1304          dst_mask_int = 1
1305        elsewhere
1306          dst_mask_int = 0
1307        endwhere
1308
1309        where (grid1_mask)
1310          src_mask_int = 1
1311        elsewhere
1312          src_mask_int = 0
1313        endwhere
1314      else
1315        allocate (src_mask_int(grid2_size), &
1316                  dst_mask_int(grid1_size))
1317
1318        where (grid1_mask)
1319          dst_mask_int = 1
1320        elsewhere
1321          dst_mask_int = 0
1322        endwhere
1323
1324        where (grid2_mask)
1325          src_mask_int = 1
1326        elsewhere
1327          src_mask_int = 0
1328        endwhere
1329      endif
1330
1331!-----------------------------------------------------------------------
1332!
1333!     change units of lat/lon coordinates if input units different
1334!     from radians. if this is the second mapping, the conversion has
1335!     alread been done.
1336!
1337!-----------------------------------------------------------------------
1338
1339      if (grid1_units(1:7) == 'degrees' .and. direction == 1) then
1340        grid1_center_lat = grid1_center_lat/deg2rad
1341        grid1_center_lon = grid1_center_lon/deg2rad
1342        grid1_corner_lat = grid1_corner_lat/deg2rad
1343        grid1_corner_lon = grid1_corner_lon/deg2rad
1344      endif
1345
1346      if (grid2_units(1:7) == 'degrees' .and. direction == 1) then
1347        grid2_center_lat = grid2_center_lat/deg2rad
1348        grid2_center_lon = grid2_center_lon/deg2rad
1349        grid2_corner_lat = grid2_corner_lat/deg2rad
1350        grid2_corner_lon = grid2_corner_lon/deg2rad
1351      endif
1352
1353!-----------------------------------------------------------------------
1354!
1355!     write mapping data
1356!
1357!-----------------------------------------------------------------------
1358
1359      if (direction == 1) then
1360        itmp1 = nc_srcgrddims_id
1361        itmp2 = nc_dstgrddims_id
1362      else
1363        itmp2 = nc_srcgrddims_id
1364        itmp1 = nc_dstgrddims_id
1365      endif
1366
1367      ncstat = nf_put_var_int(nc_file_id, itmp1, grid1_dims)
1368      call netcdf_error_handler(ncstat)
1369
1370      ncstat = nf_put_var_int(nc_file_id, itmp2, grid2_dims)
1371      call netcdf_error_handler(ncstat)
1372
1373      ncstat = nf_put_var_int(nc_file_id, nc_srcgrdimask_id,  &
1374                              src_mask_int)
1375      call netcdf_error_handler(ncstat)
1376
1377      ncstat = nf_put_var_int(nc_file_id, nc_dstgrdimask_id, &
1378                              dst_mask_int)
1379      call netcdf_error_handler(ncstat)
1380
1381      deallocate(src_mask_int, dst_mask_int)
1382
1383      if (direction == 1) then
1384        itmp1 = nc_srcgrdcntrlat_id
1385        itmp2 = nc_srcgrdcntrlon_id
1386        itmp3 = nc_srcgrdcrnrlat_id
1387        itmp4 = nc_srcgrdcrnrlon_id
1388      else
1389        itmp1 = nc_dstgrdcntrlat_id
1390        itmp2 = nc_dstgrdcntrlon_id
1391        itmp3 = nc_dstgrdcrnrlat_id
1392        itmp4 = nc_dstgrdcrnrlon_id
1393      endif
1394
1395      ncstat = nf_put_var_double(nc_file_id, itmp1, grid1_center_lat)
1396      call netcdf_error_handler(ncstat)
1397
1398      ncstat = nf_put_var_double(nc_file_id, itmp2, grid1_center_lon)
1399      call netcdf_error_handler(ncstat)
1400
1401      ncstat = nf_put_var_double(nc_file_id, itmp3, grid1_corner_lat)
1402      call netcdf_error_handler(ncstat)
1403
1404      ncstat = nf_put_var_double(nc_file_id, itmp4, grid1_corner_lon)
1405      call netcdf_error_handler(ncstat)
1406
1407      if (direction == 1) then
1408        itmp1 = nc_dstgrdcntrlat_id
1409        itmp2 = nc_dstgrdcntrlon_id
1410        itmp3 = nc_dstgrdcrnrlat_id
1411        itmp4 = nc_dstgrdcrnrlon_id
1412      else
1413        itmp1 = nc_srcgrdcntrlat_id
1414        itmp2 = nc_srcgrdcntrlon_id
1415        itmp3 = nc_srcgrdcrnrlat_id
1416        itmp4 = nc_srcgrdcrnrlon_id
1417      endif
1418
1419      ncstat = nf_put_var_double(nc_file_id, itmp1, grid2_center_lat)
1420      call netcdf_error_handler(ncstat)
1421
1422      ncstat = nf_put_var_double(nc_file_id, itmp2, grid2_center_lon)
1423      call netcdf_error_handler(ncstat)
1424
1425      ncstat = nf_put_var_double(nc_file_id, itmp3, grid2_corner_lat)
1426      call netcdf_error_handler(ncstat)
1427
1428      ncstat = nf_put_var_double(nc_file_id, itmp4, grid2_corner_lon)
1429      call netcdf_error_handler(ncstat)
1430
1431      if (direction == 1) then
1432        itmp1 = nc_srcgrdarea_id
1433        itmp2 = nc_srcgrdfrac_id
1434        itmp3 = nc_dstgrdarea_id
1435        itmp4 = nc_dstgrdfrac_id
1436      else
1437        itmp1 = nc_dstgrdarea_id
1438        itmp2 = nc_dstgrdfrac_id
1439        itmp3 = nc_srcgrdarea_id
1440        itmp4 = nc_srcgrdfrac_id
1441      endif
1442
1443      if (luse_grid1_area) then
1444        ncstat = nf_put_var_double(nc_file_id, itmp1, grid1_area_in)
1445      else
1446        ncstat = nf_put_var_double(nc_file_id, itmp1, grid1_area)
1447      endif
1448      call netcdf_error_handler(ncstat)
1449
1450      ncstat = nf_put_var_double(nc_file_id, itmp2, grid1_frac)
1451      call netcdf_error_handler(ncstat)
1452
1453      if (luse_grid2_area) then
1454        ncstat = nf_put_var_double(nc_file_id, itmp3, grid2_area)
1455      else
1456        ncstat = nf_put_var_double(nc_file_id, itmp3, grid2_area)
1457      endif
1458      call netcdf_error_handler(ncstat)
1459
1460      ncstat = nf_put_var_double(nc_file_id, itmp4, grid2_frac)
1461      call netcdf_error_handler(ncstat)
1462
1463      if (direction == 1) then
1464        ncstat = nf_put_var_int(nc_file_id, nc_srcadd_id,  &
1465                                grid1_add_map1)
1466        call netcdf_error_handler(ncstat)
1467
1468        ncstat = nf_put_var_int(nc_file_id, nc_dstadd_id,  &
1469                                grid2_add_map1)
1470        call netcdf_error_handler(ncstat)
1471
1472        if (num_wts == 1) then
1473          ncstat = nf_put_var_double(nc_file_id, nc_rmpmatrix_id,  &
1474                                     wts_map1)
1475          call netcdf_error_handler(ncstat)
1476        else
1477          allocate(wts1(num_links_map1),wts2(num_wts-1,num_links_map1))
1478
1479          wts1 = wts_map1(1,:)
1480          wts2 = wts_map1(2:,:)
1481
1482          ncstat = nf_put_var_double(nc_file_id, nc_rmpmatrix_id,  &
1483                                     wts1)
1484          call netcdf_error_handler(ncstat)
1485          ncstat = nf_put_var_double(nc_file_id, nc_rmpmatrix2_id,  &
1486                                     wts2)
1487          call netcdf_error_handler(ncstat)
1488          deallocate(wts1,wts2)
1489        endif
1490      else
1491        ncstat = nf_put_var_int(nc_file_id, nc_srcadd_id,  &
1492                                grid2_add_map2)
1493        call netcdf_error_handler(ncstat)
1494
1495        ncstat = nf_put_var_int(nc_file_id, nc_dstadd_id,  &
1496                                grid1_add_map2)
1497        call netcdf_error_handler(ncstat)
1498
1499        if (num_wts == 1) then
1500          ncstat = nf_put_var_double(nc_file_id, nc_rmpmatrix_id,  &
1501                                     wts_map2)
1502          call netcdf_error_handler(ncstat)
1503        else
1504          allocate(wts1(num_links_map2),wts2(num_wts-1,num_links_map2))
1505
1506          wts1 = wts_map2(1,:)
1507          wts2 = wts_map2(2:,:)
1508
1509          ncstat = nf_put_var_double(nc_file_id, nc_rmpmatrix_id,  &
1510                                     wts1)
1511          call netcdf_error_handler(ncstat)
1512          ncstat = nf_put_var_double(nc_file_id, nc_rmpmatrix2_id,  &
1513                                     wts2)
1514          call netcdf_error_handler(ncstat)
1515          deallocate(wts1,wts2)
1516        endif
1517      endif
1518
1519      ncstat = nf_close(nc_file_id)
1520      call netcdf_error_handler(ncstat)
1521
1522!-----------------------------------------------------------------------
1523
1524      end subroutine write_remap_csm
1525
1526!***********************************************************************
1527
1528      subroutine sort_add(add1, add2, weights)
1529
1530!-----------------------------------------------------------------------
1531!
1532!     this routine sorts address and weight arrays based on the
1533!     destination address with the source address as a secondary
1534!     sorting criterion.  the method is a standard heap sort.
1535!
1536!-----------------------------------------------------------------------
1537
1538      use kinds_mod     ! defines common data types
1539      use constants     ! defines common scalar constants
1540
1541      implicit none
1542
1543!-----------------------------------------------------------------------
1544!
1545!     Input and Output arrays
1546!
1547!-----------------------------------------------------------------------
1548
1549      integer (kind=int_kind), intent(inout), dimension(:) :: &
1550              add1,        & ! destination address array (num_links)
1551              add2        ! source      address array
1552
1553      real (kind=dbl_kind), intent(inout), dimension(:,:) :: &
1554              weights     ! remapping weights (num_wts, num_links)
1555
1556!-----------------------------------------------------------------------
1557!
1558!     local variables
1559!
1560!-----------------------------------------------------------------------
1561
1562      integer (kind=int_kind) :: &
1563                num_links,           & ! num of links for this mapping
1564                num_wts,             & ! num of weights for this mapping
1565                add1_tmp, add2_tmp,  & ! temp for addresses during swap
1566                nwgt, &
1567                lvl, final_lvl,      & ! level indexes for heap sort levels
1568                chk_lvl1, chk_lvl2, max_lvl
1569
1570      real (kind=dbl_kind), dimension(SIZE(weights,DIM=1)) :: &
1571                wgttmp              ! temp for holding wts during swap
1572
1573!-----------------------------------------------------------------------
1574!
1575!     determine total number of links to sort and number of weights
1576!
1577!-----------------------------------------------------------------------
1578
1579      num_links = SIZE(add1)
1580      num_wts   = SIZE(weights, DIM=1)
1581
1582!-----------------------------------------------------------------------
1583!
1584!     start at the lowest level (N/2) of the tree and sift lower
1585!     values to the bottom of the tree, promoting the larger numbers
1586!
1587!-----------------------------------------------------------------------
1588
1589      do lvl=num_links/2,1,-1
1590
1591        final_lvl = lvl
1592        add1_tmp = add1(lvl)
1593        add2_tmp = add2(lvl)
1594        wgttmp(:) = weights(:,lvl)
1595
1596        !***
1597        !*** loop until proper level is found for this link, or reach
1598        !*** bottom
1599        !***
1600
1601        sift_loop1: do
1602
1603          !***
1604          !*** find the largest of the two daughters
1605          !***
1606
1607          chk_lvl1 = 2*final_lvl
1608          chk_lvl2 = 2*final_lvl+1
1609          if (chk_lvl1 .EQ. num_links) chk_lvl2 = chk_lvl1
1610
1611          if ((add1(chk_lvl1) >  add1(chk_lvl2)) .OR. &
1612             ((add1(chk_lvl1) == add1(chk_lvl2)) .AND. &
1613              (add2(chk_lvl1) >  add2(chk_lvl2)))) then
1614            max_lvl = chk_lvl1
1615          else
1616            max_lvl = chk_lvl2
1617          endif
1618
1619          !***
1620          !*** if the parent is greater than both daughters,
1621          !*** the correct level has been found
1622          !***
1623
1624          if ((add1_tmp .GT. add1(max_lvl)) .OR. &
1625             ((add1_tmp .EQ. add1(max_lvl)) .AND. &
1626              (add2_tmp .GT. add2(max_lvl)))) then
1627            add1(final_lvl) = add1_tmp
1628            add2(final_lvl) = add2_tmp
1629            weights(:,final_lvl) = wgttmp(:)
1630            exit sift_loop1
1631
1632          !***
1633          !*** otherwise, promote the largest daughter and push
1634          !*** down one level in the tree.  if haven't reached
1635          !*** the end of the tree, repeat the process.  otherwise
1636          !*** store last values and exit the loop
1637          !***
1638
1639          else
1640            add1(final_lvl) = add1(max_lvl)
1641            add2(final_lvl) = add2(max_lvl)
1642            weights(:,final_lvl) = weights(:,max_lvl)
1643
1644            final_lvl = max_lvl
1645            if (2*final_lvl > num_links) then
1646              add1(final_lvl) = add1_tmp
1647              add2(final_lvl) = add2_tmp
1648              weights(:,final_lvl) = wgttmp(:)
1649              exit sift_loop1
1650            endif
1651          endif
1652        end do sift_loop1
1653      end do
1654
1655!-----------------------------------------------------------------------
1656!
1657!     now that the heap has been sorted, strip off the top (largest)
1658!     value and promote the values below
1659!
1660!-----------------------------------------------------------------------
1661
1662      do lvl=num_links,3,-1
1663
1664        !***
1665        !*** move the top value and insert it into the correct place
1666        !***
1667
1668        add1_tmp = add1(lvl)
1669        add1(lvl) = add1(1)
1670
1671        add2_tmp = add2(lvl)
1672        add2(lvl) = add2(1)
1673
1674        wgttmp(:) = weights(:,lvl)
1675        weights(:,lvl) = weights(:,1)
1676
1677        !***
1678        !*** as above this loop sifts the tmp values down until proper
1679        !*** level is reached
1680        !***
1681
1682        final_lvl = 1
1683
1684        sift_loop2: do
1685
1686          !***
1687          !*** find the largest of the two daughters
1688          !***
1689
1690          chk_lvl1 = 2*final_lvl
1691          chk_lvl2 = 2*final_lvl+1
1692          if (chk_lvl2 >= lvl) chk_lvl2 = chk_lvl1
1693
1694          if ((add1(chk_lvl1) >  add1(chk_lvl2)) .OR. &
1695             ((add1(chk_lvl1) == add1(chk_lvl2)) .AND. &
1696              (add2(chk_lvl1) >  add2(chk_lvl2)))) then
1697            max_lvl = chk_lvl1
1698          else
1699            max_lvl = chk_lvl2
1700          endif
1701
1702          !***
1703          !*** if the parent is greater than both daughters,
1704          !*** the correct level has been found
1705          !***
1706
1707          if ((add1_tmp >  add1(max_lvl)) .OR. &
1708             ((add1_tmp == add1(max_lvl)) .AND. &
1709              (add2_tmp >  add2(max_lvl)))) then
1710            add1(final_lvl) = add1_tmp
1711            add2(final_lvl) = add2_tmp
1712            weights(:,final_lvl) = wgttmp(:)
1713            exit sift_loop2
1714
1715          !***
1716          !*** otherwise, promote the largest daughter and push
1717          !*** down one level in the tree.  if haven't reached
1718          !*** the end of the tree, repeat the process.  otherwise
1719          !*** store last values and exit the loop
1720          !***
1721
1722          else
1723            add1(final_lvl) = add1(max_lvl)
1724            add2(final_lvl) = add2(max_lvl)
1725            weights(:,final_lvl) = weights(:,max_lvl)
1726
1727            final_lvl = max_lvl
1728            if (2*final_lvl >= lvl) then
1729              add1(final_lvl) = add1_tmp
1730              add2(final_lvl) = add2_tmp
1731              weights(:,final_lvl) = wgttmp(:)
1732              exit sift_loop2
1733            endif
1734          endif
1735        end do sift_loop2
1736      end do
1737
1738      !***
1739      !*** swap the last two entries
1740      !***
1741
1742
1743      add1_tmp = add1(2)
1744      add1(2)  = add1(1)
1745      add1(1)  = add1_tmp
1746
1747      add2_tmp = add2(2)
1748      add2(2)  = add2(1)
1749      add2(1)  = add2_tmp
1750
1751      wgttmp (:)   = weights(:,2)
1752      weights(:,2) = weights(:,1)
1753      weights(:,1) = wgttmp (:)
1754
1755!-----------------------------------------------------------------------
1756
1757      end subroutine sort_add
1758
1759!***********************************************************************
1760
1761      end module remap_write
1762
1763!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.