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/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/WEIGHTS/src – NEMO

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/WEIGHTS/src/scripinterp_mod.F90 @ 9295

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

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

File size: 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.