New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
scripinterp_mod.F90 in branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/TOOLS/WEIGHTS/src – NEMO

source: branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/TOOLS/WEIGHTS/src/scripinterp_mod.F90 @ 5312

Last change on this file since 5312 was 5312, checked in by timgraham, 9 years ago

Reset svn:keywords Id property

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