source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/TOOLS/WEIGHTS/nocsutil/scripgrid_mod.F90 @ 5550

Last change on this file since 5550 was 2718, checked in by sga, 10 years ago

NEMO trunk: small corrections to the weights generation code in the TOOLS directory for nemo 3.3

  • typos in README file corrected
  • scripgrid.F90 now allows any input nemo grid filename (but lat/lon names have to begin with glam/gphi)
  • scripinterp.F90 now actually does something
File size: 28.2 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     This module creates grid description files for input to the SCRIP code
4!
5!-----------------------------------------------------------------------
6
7MODULE scripgrid_mod
8
9  USE kinds_mod
10  USE constants
11  USE iounits
12  USE netcdf
13  USE netcdf_mod
14 
15  IMPLICIT NONE
16
17  !-----------------------------------------------------------------------
18  !     module variables that describe the grid
19
20  INTEGER (kind=int_kind), parameter :: &
21       grid_rank = 2,         &
22       grid_corners = 4
23  INTEGER (kind=int_kind) :: nx, ny, grid_size
24  INTEGER (kind=int_kind), dimension(2) :: &
25       grid_dims,    &  ! size of x, y dimensions
26       grid_dim_ids     ! ids of the x, y dimensions
27  INTEGER (kind=int_kind), ALLOCATABLE, DIMENSION(:) :: &
28       grid_imask          ! land-sea mask
29  REAL (kind=int_kind), ALLOCATABLE, DIMENSION(:) :: &
30       grid_center_lat, &  ! lat/lon coordinates for
31       grid_center_lon     ! each grid center in degrees
32  REAL (kind=dbl_kind), ALLOCATABLE, DIMENSION(:,:) :: &
33       grid_corner_lat, &  ! lat/lon coordinates for
34       grid_corner_lon     ! each grid corner in degrees
35  REAL (kind=dbl_kind), ALLOCATABLE, DIMENSION(:,:,:) :: &
36       corner_lon, &
37       corner_lat
38  REAL (kind=dbl_kind), PARAMETER :: circle = 360.0
39
40  !-----------------------------------------------------------------------
41  !     module variables that describe the netcdf file
42
43  INTEGER (kind=int_kind) :: &
44       ncstat,           &   ! general netCDF status variable
45       ncid_in
46
47CONTAINS
48
49  ! ==============================================================================
50 
51  SUBROUTINE convert(nm_in)
52 
53    ! -----------------------------------------------------------------------------
54    ! - input variables
55 
56    CHARACTER(char_len), INTENT(in) ::  &
57      nm_in
58 
59    ! -----------------------------------------------------------------------------
60    ! - local variables
61 
62    CHARACTER(char_len) ::  &
63      nemo_file, input_file, method, input_lon, input_lat, datagrid_file, &
64      nemogrid_file, nemo_lon, nemo_lat, corn_lon, corn_lat, nemo_mask, input_mask
65    INTEGER (kind=int_kind), dimension(2) :: &
66      offset
67    INTEGER (kind=int_kind) :: &
68      iunit, nemo_mask_value, input_mask_value
69 
70    namelist /grid_inputs/ nemo_file, input_file, datagrid_file, nemogrid_file,  &
71                           method, input_lon, input_lat, nemo_lon, nemo_lat,     &
72                           nemo_mask, nemo_mask_value, input_mask, input_mask_value
73 
74    !-----------------------------------------------------------------------
75    ! - namelist describing the processing
76    !   note that mask_value is the minimum good value,
77    !   so that where the mask is less than the value is masked
78
79    nemo_file = "coordinates.nc"
80    nemo_lon = "glamt"
81    nemo_lat = "gphit"
82    input_lon = "lon"
83    input_lat = "lat"
84    input_mask = "none"
85    input_mask_value = 0
86    datagrid_file = 'remap_data_grid.nc'
87    nemogrid_file = 'remap_nemo_grid.nc'
88
89    call get_unit(iunit)
90    open(iunit, file=nm_in, status='old', form='formatted')
91    read(iunit, nml=grid_inputs) 
92    call release_unit(iunit)
93
94    if (nemo_lon(1:4) .ne. 'glam' .or. nemo_lat(1:4) .ne. 'gphi') then
95      write(6,*) 'lon name does not start with "glam" or lat name does not start with "gphi"'
96      stop
97    endif
98
99    ! set up the names of the corner variables for a given input
100    ! the offset represents what needs to be added to (i,j) to get to the correct
101    ! element in the corner arrays to correspond to the point northeast of the center
102    if (nemo_lon(5:5) == "t") then
103      corn_lon = "glamf"
104      corn_lat = "gphif"
105      offset = (/ 0,0 /)
106    else if (nemo_lon(5:5) == "u") then
107      corn_lon = "glamv"
108      corn_lat = "gphiv"
109      offset = (/ 1,0 /)
110    else if (nemo_lon(5:5) == "v") then
111      corn_lon = "glamu"
112      corn_lat = "gphiu"
113      offset = (/ 0,1 /)
114    else
115      write(6,*) 'unknown nemo_lon name'
116      stop
117    endif
118
119    write(6,*) "processing " // trim(nemo_file)
120    call convertNEMO(nemo_file, nemo_lon, nemo_lat, corn_lon, corn_lat, &
121                     offset, nemogrid_file)
122
123    write(6,*) "processing regular grid"
124    call convertFLUX(input_file, input_lon, input_lat, &
125                     input_mask, input_mask_value, datagrid_file)
126
127  END SUBROUTINE convert
128 
129  ! ==============================================================================
130 
131  SUBROUTINE convertNEMO(grid_file_in, cent_lon, cent_lat, corn_lon, corn_lat, &
132                         off, grid_file_out)
133 
134    !-----------------------------------------------------------------------
135    !
136    !     This routine converts a NEMO coordinates.nc file to a remapping grid file.
137    !
138   
139    CHARACTER(char_len), INTENT(in) :: cent_lon, cent_lat, corn_lon, corn_lat
140    INTEGER (kind=int_kind), INTENT(in), DIMENSION(2) :: off
141    CHARACTER(char_len), INTENT(in) :: grid_file_out
142    CHARACTER(char_len), INTENT(in) :: grid_file_in
143
144    !-----------------------------------------------------------------------
145    !     module variables that describe the grid
146 
147    CHARACTER(char_len), parameter ::  &
148         grid_name = 'Remapped NEMO grid for SCRIP'
149 
150    !-----------------------------------------------------------------------
151    !     grid coordinates and masks
152 
153    REAL (kind=dbl_kind), ALLOCATABLE, DIMENSION(:,:) :: &
154         clon, clat, &             ! expanded corner arrays
155         glam, &                   ! center longitude
156         gphi, &                   ! center latitude
157         glamc, &                  ! corner longitude
158         gphic                     ! corner latitude
159 
160    !-----------------------------------------------------------------------
161    !     other local variables
162 
163    INTEGER (kind=int_kind) :: i, j, n, iunit, im1, jm1, imid, isame, ic, jc
164    INTEGER (kind=int_kind) :: varid_lam, varid_phi, varid_lamc, varid_phic
165    INTEGER (kind=int_kind) :: jdim
166    INTEGER (kind=int_kind), dimension(4) :: grid_dimids  ! input fields have 4 dims
167    REAL (kind=dbl_kind) :: tmplon, dxt, dyt
168 
169    !-----------------------------------------------------------------------
170    !     read in grid info
171    !
172    !     For NEMO input grids, assume that variable names are glam, glamc etc.
173    !     Assume that 1st 2 dimensions of these variables are x and y directions.
174    !     These assumptions are made by NEMO, so should be valid for coordinates.nc.
175    !
176    ! write in nf90 calls (without error handling) and then think about
177    ! making more readable by taking chunks into ncutil
178    !
179 
180    ncstat = nf90_open( grid_file_in, NF90_NOWRITE, ncid_in )
181    call netcdf_error_handler(ncstat)
182 
183    ! find dimids for 'glam'
184    ! use dimids to get dimlengths
185    ! allocate glam array
186    ! get glam from file
187 
188    ncstat = nf90_inq_varid( ncid_in, cent_lon, varid_lam )
189    call netcdf_error_handler(ncstat)
190    ncstat = nf90_inq_varid( ncid_in, corn_lon, varid_lamc )
191    call netcdf_error_handler(ncstat)
192    ncstat = nf90_inq_varid( ncid_in, cent_lat, varid_phi )
193    call netcdf_error_handler(ncstat)
194    ncstat = nf90_inq_varid( ncid_in, corn_lat, varid_phic )
195    call netcdf_error_handler(ncstat)
196 
197    ncstat = nf90_inquire_variable( ncid_in, varid_lam, dimids=grid_dimids(:) )
198    call netcdf_error_handler(ncstat)
199    DO jdim = 1, SIZE(grid_dims)
200       ncstat = nf90_inquire_dimension( ncid_in, grid_dimids(jdim),  &
201                                        len=grid_dims(jdim) )
202       call netcdf_error_handler(ncstat)
203    END DO
204    nx = grid_dims(1)
205    ny = grid_dims(2)
206    grid_size = nx * ny
207    WRITE(*,FMT='("Input grid dimensions are:",2i6)') nx, ny
208   
209    ! assume that dimensions are all the same as glam
210    ALLOCATE( glam(nx,ny), glamc(nx,ny), gphi(nx,ny), gphic(nx,ny) )
211    ncstat = nf90_get_var( ncid_in, varid_lam, glam(:,:) )
212    call netcdf_error_handler(ncstat)
213    ncstat = nf90_get_var( ncid_in, varid_lamc, glamc(:,:) )
214    call netcdf_error_handler(ncstat)
215    ncstat = nf90_get_var( ncid_in, varid_phi, gphi(:,:) )
216    call netcdf_error_handler(ncstat)
217    ncstat = nf90_get_var( ncid_in, varid_phic, gphic(:,:) )
218    call netcdf_error_handler(ncstat)
219   
220    !-----------------------------------------------------------------------
221    ! - Mask is all ocean for now
222 
223    ALLOCATE( grid_imask(grid_size) )
224    grid_imask(:) = 1
225 
226    !-----------------------------------------------------------------------
227    ! corners are arranged as follows:    4 3
228    !                                     1 2
229    !
230    ! Assume that cyclic grids have 2 wrap columns in coordinates.nc
231    !  (this is the case for ORCA grids)
232    !
233 
234    ! -----------------------------------------------------------------------------
235    ! create a single pair of arrays for the corners where clon(1,1) corresponds
236    ! to the south west corner of a box containing glam(1,1)
237    ! various special cases then apply
238    !   bottom row:  assume clon(:,j) = clon(:,j+1)
239 
240    ALLOCATE ( clon(nx+1,ny+1), clat(nx+1,ny+1) )
241
242    ! first the easy internal points
243    DO j = 2,ny
244      DO i = 2,nx
245        ic = i + off(1) - 1
246        jc = j + off(2) - 1
247        clon(i,j) = glamc(ic,jc)
248        clat(i,j) = gphic(ic,jc)
249      ENDDO
250    ENDDO
251
252    ! then the tricky boundary points
253    imid = (nx-1)/2 + 1
254    DO j = 1,ny+1,ny
255      DO i = 1,nx+1,nx
256        ic = i + off(1) - 1
257        jc = j + off(2) - 1
258        if (ic == 0 .and. jc == 0) then
259          clon(i,j) = glamc(nx,1)
260          clat(i,j) = gphic(nx,1) - (gphic(nx,2)-gphic(nx,1))
261        else if (ic == nx+1 .and. jc == 0) then
262          clon(i,j) = glamc(1,1)
263          clat(i,j) = gphic(1,1) - (gphic(1,2)-gphic(1,1))
264        else if (ic == 0 .and. jc == ny+1) then
265          isame = 2*imid - nx + 1
266          clon(i,j) = glamc(isame,jc-1)
267          clat(i,j) = gphic(isame,jc-1)
268        else if (ic == nx+1 .and. jc == ny+1) then
269          isame = 2*imid
270          clon(i,j) = glamc(isame,jc-1)
271          clat(i,j) = gphic(isame,jc-1)
272        else if (ic == 0) then
273          clon(i,j) = glamc(nx,jc)
274          clat(i,j) = gphic(nx,jc)
275        else if (jc == 0) then
276          clon(i,j) = glamc(ic,1)
277          clat(i,j) = gphic(ic,1) - (gphic(ic,2)-gphic(ic,1))
278        else if (ic == nx+1) then
279          clon(i,j) = glamc(1,jc)
280          clat(i,j) = gphic(1,jc)
281        else if (jc == ny+1) then
282          isame = 2*imid - ic + 1
283          clon(i,j) = glamc(isame,jc-1)
284          clat(i,j) = gphic(isame,jc-1)
285        endif
286      ENDDO
287    ENDDO
288
289    ALLOCATE ( corner_lon(4,nx,ny), corner_lat(4,nx,ny) )
290 
291    ! top-right corner
292    corner_lon(3,:,:) = clon(2:nx+1,2:ny+1)
293    corner_lat(3,:,:) = clat(2:nx+1,2:ny+1)
294 
295    ! top-left corner
296    corner_lon(4,:,:) = clon(1:nx,2:ny+1)
297    corner_lat(4,:,:) = clat(1:nx,2:ny+1)
298 
299    ! bottom-right corner
300    corner_lon(2,:,:) = clon(2:nx+1,1:ny)
301    corner_lat(2,:,:) = clat(2:nx+1,1:ny)
302 
303    ! bottom-left corner
304    corner_lon(1,:,:) = clon(1:nx,1:ny)
305    corner_lat(1,:,:) = clat(1:nx,1:ny)
306 
307  ! For [N, E, W]-ward extrapolation near the poles, should we use stereographic (or
308  ! similar) projection?  This issue will come for V,F interpolation, and for all
309  ! grids with non-cyclic grids.
310 
311    ! -----------------------------------------------------------------------------
312    !     correct for 0,2pi longitude crossings
313    !      (In practice this means putting all corners into 0,2pi range
314    !       and ensuring that no box corners are miles from each other.
315    !       3pi/2 is used as threshold - I think this is quite arbitrary.)
316 
317    corner_lon(:,:,:) = MODULO( corner_lon(:,:,:), circle )
318    DO n = 2, grid_corners
319       WHERE    ( corner_lon(n,:,:) - corner_lon(n-1,:,:) < -three*circle*0.25 )
320          corner_lon(n,:,:) = corner_lon(n,:,:) + circle
321       ELSEWHERE( corner_lon(n,:,:) - corner_lon(n-1,:,:) >  three*circle*0.25 )
322          corner_lon(n,:,:) = corner_lon(n,:,:) - circle
323       END WHERE
324    END DO
325 
326    ! -----------------------------------------------------------------------------
327    ! - put longitudes on smooth grid
328 
329 !  call mouldlon(glam,nx,ny)
330 !  call mouldlon(corner_lon(1,:,:),nx,ny)
331 !  call mouldlon(corner_lon(2,:,:),nx,ny)
332 !  call mouldlon(corner_lon(3,:,:),nx,ny)
333 !  call mouldlon(corner_lon(4,:,:),nx,ny)
334
335    ! -----------------------------------------------------------------------------
336    ! - reshape for SCRIP input format
337 
338    ALLOCATE( grid_center_lon(grid_size), grid_center_lat(grid_size) )
339   
340    grid_center_lon(:) = RESHAPE( glam(:,:), (/ grid_size /) )
341    grid_center_lat(:) = RESHAPE( gphi(:,:), (/ grid_size /) )
342 
343    DEALLOCATE( glam, gphi, glamc, gphic )
344 
345    ALLOCATE( grid_corner_lon(4, grid_size), grid_corner_lat(4, grid_size) )
346 
347    grid_corner_lon(:,:) = RESHAPE( corner_lon(:,:,:), (/ 4, grid_size /) )
348    grid_corner_lat(:,:) = RESHAPE( corner_lat(:,:,:), (/ 4, grid_size /) )
349 
350    DEALLOCATE( corner_lon, corner_lat )
351 
352    CALL createSCRIPgrid(grid_file_out, grid_name)
353 
354  END SUBROUTINE convertNEMO
355 
356  ! ==============================================================================
357 
358  SUBROUTINE convertFLUX(grid_file_in, name_lon, name_lat, &
359                         name_mask, value_mask, grid_file_out)
360 
361    !-----------------------------------------------------------------------
362    !
363    !     This routine creates a remapping grid file from an input grid.
364    !
365    !-----------------------------------------------------------------------
366   
367    CHARACTER(char_len), INTENT(in) ::  &
368         grid_file_in, name_lon, name_lat, name_mask, grid_file_out
369    INTEGER (kind=int_kind) :: value_mask
370 
371    !-----------------------------------------------------------------------
372    !     variables that describe the grid
373 
374    CHARACTER(char_len), parameter ::  &
375         grid_name = 'Remapped regular grid for SCRIP'
376 
377    !-----------------------------------------------------------------------
378    !     grid coordinates (note that a flux file just has lon and lat)
379 
380    REAL (kind=dbl_kind), ALLOCATABLE, DIMENSION(:) :: &
381         lam, phi
382    REAL (kind=dbl_kind), ALLOCATABLE, DIMENSION(:,:) :: &
383         glam, &                   ! longitude
384         gphi, &                   ! latitude
385         glamc, &
386         gphic
387    REAL (kind=dbl_kind), ALLOCATABLE, DIMENSION(:,:) :: mask
388 
389    !-----------------------------------------------------------------------
390    !     other local variables
391 
392    INTEGER (kind=int_kind) :: i, j, n, iunit, im1, jm1
393    INTEGER (kind=int_kind) :: varid_lam, varid_phi, varid_mask
394    INTEGER (kind=int_kind) :: jdim, nspace
395    INTEGER (kind=int_kind), dimension(4) :: grid_dimids  ! input fields have 4 dims
396    REAL (kind=dbl_kind) :: tmplon, dxt, dyt
397 
398    !-----------------------------------------------------------------------
399    !     read in grid info
400    !
401    !     For NEMO input grids, assume that variable names are glam, glamc etc.
402    !     Assume that 1st 2 dimensions of these variables are x and y directions.
403    !     These assumptions are made by NEMO, so should be valid for coordinates.nc.
404    !
405    ! write in nf90 calls (without error handling) and then think about
406    ! making more readable by taking chunks into ncutil
407 
408    ncstat = nf90_open( grid_file_in, NF90_NOWRITE, ncid_in )
409    call netcdf_error_handler(ncstat)
410 
411    ! find dimids for 'glamt'
412    ! use dimids to get dimlengths
413    ! allocate glam array
414    ! get glam from file
415 
416    ncstat = nf90_inq_varid( ncid_in, name_lat, varid_phi )
417    call netcdf_error_handler(ncstat)
418    ncstat = nf90_inq_varid( ncid_in, name_lon, varid_lam )
419    call netcdf_error_handler(ncstat)
420 
421    ncstat = nf90_inquire_variable( ncid_in, varid_lam, ndims=nspace )
422    call netcdf_error_handler(ncstat)
423
424    if (nspace == 1) then
425      ncstat = nf90_inquire_variable( ncid_in, varid_lam, dimids=grid_dimids(:1) )
426      call netcdf_error_handler(ncstat)
427      ncstat = nf90_inquire_variable( ncid_in, varid_phi, dimids=grid_dimids(2:) )
428      call netcdf_error_handler(ncstat)
429      ncstat = nf90_inquire_dimension( ncid_in, grid_dimids(1), len=grid_dims(1) )
430      call netcdf_error_handler(ncstat)
431      ncstat = nf90_inquire_dimension( ncid_in, grid_dimids(2), len=grid_dims(2) )
432      call netcdf_error_handler(ncstat)
433      nx = grid_dims(1)
434      ny = grid_dims(2)
435      grid_size = nx * ny
436      WRITE(*,FMT='("Input grid dimensions are:",2i6)') nx, ny
437   
438      ALLOCATE( lam(nx), phi(ny) )
439      write(6,*) 'double'
440      ncstat = nf90_get_var( ncid_in, varid_lam, lam )
441      call netcdf_error_handler(ncstat)
442      ncstat = nf90_get_var( ncid_in, varid_phi, phi )
443      call netcdf_error_handler(ncstat)
444   
445      ALLOCATE( glam(nx,ny), gphi(nx,ny))
446      write(6,*) shape(lam),shape(phi)
447      glam(:,:) = SPREAD(lam,2,ny)
448      gphi(:,:) = SPREAD(phi,1,nx)
449    else
450
451      ncstat = nf90_inquire_variable( ncid_in, varid_lam, dimids=grid_dimids(:2) )
452      call netcdf_error_handler(ncstat)
453      ncstat = nf90_inquire_dimension( ncid_in, grid_dimids(1), len=grid_dims(1) )
454      call netcdf_error_handler(ncstat)
455      ncstat = nf90_inquire_dimension( ncid_in, grid_dimids(2), len=grid_dims(2) )
456      call netcdf_error_handler(ncstat)
457      nx = grid_dims(1)
458      ny = grid_dims(2)
459      grid_size = nx * ny
460      WRITE(*,FMT='("Input grid dimensions are:",2i6)') nx, ny
461   
462      ALLOCATE( glam(nx,ny), gphi(nx,ny))
463      ncstat = nf90_get_var( ncid_in, varid_lam, glam )
464      call netcdf_error_handler(ncstat)
465      ncstat = nf90_get_var( ncid_in, varid_phi, gphi )
466      call netcdf_error_handler(ncstat)
467
468    endif
469    write(6,*) grid_size,nx,ny
470 
471    ALLOCATE(glamc(0:nx,0:ny), gphic(0:nx,0:ny) )
472
473    ! - for now a simple average to get top right box corners
474    ! - glamc(i,j), gphic(i,j) are top right coordinates of box containing
475    ! - glam(i,j),gphi(i,j)
476    write(6,*) 'averaging'
477    write(6,*) size(gphic),size(gphi)
478    gphic(1:nx,1:ny-1) = 0.5*(gphi(:,1:ny-1)+gphi(:,2:ny))
479    write(6,*) size(glamc),size(glam)
480    glamc(1:nx-1,1:ny) = 0.5*(glam(1:nx-1,:)+glam(2:nx,:))
481 
482    ! - left and right column of longitudes
483    write(6,*) 'columns'
484    glamc(nx,1:ny) = 1.5*glam(nx,:)-0.5*glam(nx-1,:)
485    glamc( 0,1:ny) = 1.5*glam(1,:)-0.5*glam(2,:)
486    glamc(nx, 0) = glamc(nx,1)
487    glamc( 0, 0) = glamc( 0,1)
488
489    ! - top and bottom row of latitudes by extrapolation
490    write(6,*) 'rows'
491    gphic(1:nx,ny) = 1.5*gphi(:,ny)-0.5*gphi(:,ny-1)
492    gphic(1:nx, 0) = 1.5*gphi(:,1)-0.5*gphi(:,2)
493    gphic( 0,ny) = gphic(1,ny)
494    gphic( 0, 0) = gphic(1, 0)
495
496    !-----------------------------------------------------------------------
497 
498    write(6,*) 'allocating'
499    ALLOCATE( grid_imask(grid_size) )
500    grid_imask(:) = 1
501    write(6,*) name_mask
502    if (trim(name_mask) /= "none") then
503      write(6,*) 'masking'
504      ncstat = nf90_inq_varid( ncid_in, name_mask, varid_mask )
505      call netcdf_error_handler(ncstat)
506      ALLOCATE( mask(nx,ny) )
507      write(6,*) 'reading mask'
508      ncstat = nf90_get_var( ncid_in, varid_mask, mask )
509      call netcdf_error_handler(ncstat)
510      write(6,*) 'setting mask'
511      WHERE ( RESHAPE(mask(:,:),(/ grid_size /)) < value_mask)
512        grid_imask = 0
513      END WHERE
514      write(6,*) 'masked'
515    END IF
516 
517    !-----------------------------------------------------------------------
518    ! corners are arranged as follows:    4 3
519    !                                     1 2
520 
521    ALLOCATE ( corner_lon(4,nx,ny), corner_lat(4,nx,ny) )
522 
523    ! - bottom-left corner
524    corner_lon(1,:,:) = glamc(0:nx-1, 0:ny-1 )
525    corner_lat(1,:,:) = gphic(0:nx-1, 0:ny-1 )
526 
527    ! - bottom-right corner
528    corner_lon(2,:,:) = glamc(1:nx, 0:ny-1 )
529    corner_lat(2,:,:) = gphic(1:nx, 0:ny-1 )
530 
531    ! - top-right corner
532    corner_lon(3,:,:) = glamc(1:nx,1:ny)
533    corner_lat(3,:,:) = gphic(1:nx,1:ny)
534    write(6,*) corner_lat(3,nx-2:nx,ny)
535 
536    ! - top-left corner
537    corner_lon(4,:,:) = glamc(0:nx-1, 1:ny )
538    corner_lat(4,:,:) = gphic(0:nx-1, 1:ny )
539 
540  ! For [N, E, W]-ward extrapolation near the poles, should we use stereographic (or
541  ! similar) projection?  This issue will come for V,F interpolation, and for all
542  ! grids with non-cyclic grids.
543 
544    ! -----------------------------------------------------------------------------
545    !     correct for 0,2pi longitude crossings
546    !      (In practice this means putting all corners into 0,2pi range
547    !       and ensuring that no box corners are miles from each other.
548    !       3pi/2 is used as threshold - I think this is quite arbitrary.)
549 
550  ! corner_lon(:,:,:) = MODULO( corner_lon(:,:,:), circle )
551  ! DO n = 2, grid_corners
552  !    WHERE    ( corner_lon(n,:,:) - corner_lon(n-1,:,:) < -three*circle*0.25 )
553  !       corner_lon(n,:,:) = corner_lon(n,:,:) + circle
554  !    ELSEWHERE( corner_lon(n,:,:) - corner_lon(n-1,:,:) >  three*circle*0.25 )
555  !       corner_lon(n,:,:) = corner_lon(n,:,:) - circle
556  !    END WHERE
557  ! END DO
558 
559    ! -----------------------------------------------------------------------------
560    ! - reshape for SCRIP input format
561 
562    ALLOCATE( grid_center_lon(grid_size), grid_center_lat(grid_size) )
563   
564    grid_center_lon(:) = RESHAPE( glam(:,:), (/ grid_size /) )
565    grid_center_lat(:) = RESHAPE( gphi(:,:), (/ grid_size /) )
566 
567    DEALLOCATE( glam, gphi, glamc, gphic )
568 
569    ALLOCATE( grid_corner_lon(4, grid_size), grid_corner_lat(4, grid_size) )
570 
571    grid_corner_lon(:,:) = RESHAPE( corner_lon(:,:,:), (/ 4, grid_size /) )
572    grid_corner_lat(:,:) = RESHAPE( corner_lat(:,:,:), (/ 4, grid_size /) )
573 
574    DEALLOCATE( corner_lon, corner_lat )
575 
576    CALL createSCRIPgrid(grid_file_out, grid_name)
577 
578  END SUBROUTINE convertFLUX
579 
580  ! ==============================================================================
581 
582  SUBROUTINE mouldlon(lon_grid, nx, ny)
583
584    ! -----------------------------------------------------------------------------
585    ! - input variables
586
587    INTEGER, INTENT(in) :: nx, ny
588    REAL (kind=dbl_kind), INTENT(inout), DIMENSION(nx,ny) ::  &
589      lon_grid
590 
591    ! -----------------------------------------------------------------------------
592    ! - local variables
593
594    INTEGER :: ix, iy
595    REAL (kind=dbl_kind), DIMENSION(:,:), ALLOCATABLE ::  &
596      dlon
597    REAL :: step
598
599    ! -----------------------------------------------------------------------------
600    ! - try to eliminate any 360 degree steps in a grid of longitudes
601
602    ALLOCATE(dlon(nx,ny))
603
604    step = 0.75*circle
605    dlon(:,:) = 0
606    dlon(2:,:) = lon_grid(2:,:) - lon_grid(:nx-1,:)
607    WHERE (dlon > -step .AND. dlon < step)
608      dlon = 0.0
609    ELSEWHERE
610      dlon = -SIGN(circle,dlon)
611    END WHERE
612
613    ! - close your eyes this is nasty
614    DO ix = 2,nx
615      dlon(ix,:) = dlon(ix,:) + dlon(ix-1,:)
616    END DO
617    lon_grid = lon_grid + dlon
618
619  END SUBROUTINE mouldlon
620 
621  ! ==============================================================================
622 
623  SUBROUTINE createSCRIPgrid(grid_file_out, grid_name)
624 
625    ! -----------------------------------------------------------------------------
626    ! - input variables
627 
628    CHARACTER(char_len), INTENT(in) ::  &
629         grid_name, grid_file_out
630 
631    ! -----------------------------------------------------------------------------
632    ! - local variables that describe the netcdf file
633 
634    INTEGER (kind=int_kind) :: &
635         nc_grid_id,       &   ! netCDF grid dataset id
636         nc_gridsize_id,   &   ! netCDF grid size dim id
637         nc_gridcorn_id,   &   ! netCDF grid corner dim id
638         nc_gridrank_id,   &   ! netCDF grid rank dim id
639         nc_griddims_id,   &   ! netCDF grid dimensions id
640         nc_grdcntrlat_id, &   ! netCDF grid center lat id
641         nc_grdcntrlon_id, &   ! netCDF grid center lon id
642         nc_grdimask_id,   &   ! netCDF grid mask id
643         nc_gridarea_id,   &   ! netCDF grid area id
644         nc_grdcrnrlat_id, &   ! netCDF grid corner lat id
645         nc_grdcrnrlon_id      ! netCDF grid corner lon id
646 
647    ! -----------------------------------------------------------------------------
648    ! - create netCDF dataset for this grid
649    ! - rewrite in nf90
650    ! - (bring out functional blocks into ncclear for readability)
651 
652    ncstat = nf90_create (grid_file_out, NF90_CLOBBER, nc_grid_id)
653    call netcdf_error_handler(ncstat)
654    ncstat = nf90_put_att (nc_grid_id, NF90_GLOBAL, 'title', grid_name)
655    call netcdf_error_handler(ncstat)
656 
657    ! -----------------------------------------------------------------------------
658    ! - define grid size dimension
659 
660    ncstat = nf90_def_dim (nc_grid_id, 'grid_size', grid_size, nc_gridsize_id)
661    call netcdf_error_handler(ncstat)
662 
663    ! -----------------------------------------------------------------------------
664    ! - define grid rank dimension
665 
666    ncstat = nf90_def_dim (nc_grid_id, 'grid_rank', grid_rank, nc_gridrank_id)
667    call netcdf_error_handler(ncstat)
668 
669    ! -----------------------------------------------------------------------------
670    ! - define grid corner dimension
671 
672    ncstat = nf90_def_dim (nc_grid_id, 'grid_corners', grid_corners, nc_gridcorn_id)
673    call netcdf_error_handler(ncstat)
674 
675    ! -----------------------------------------------------------------------------
676    ! - define grid dim size array
677 
678    ncstat = nf90_def_var(nc_grid_id, 'grid_dims', NF90_INT, nc_gridrank_id, nc_griddims_id)
679    call netcdf_error_handler(ncstat)
680 
681    ! -----------------------------------------------------------------------------
682    ! - define grid mask
683 
684    ncstat = nf90_def_var(nc_grid_id, 'grid_imask', NF90_INT, &
685                          nc_gridsize_id, nc_grdimask_id)
686    call netcdf_error_handler(ncstat)
687    ncstat = nf90_put_att(nc_grid_id, nc_grdimask_id, 'units', 'unitless')
688    call netcdf_error_handler(ncstat)
689 
690    ! -----------------------------------------------------------------------------
691    ! - define grid center latitude array
692 
693    ncstat = nf90_def_var(nc_grid_id, 'grid_center_lat', NF90_DOUBLE, &
694                          nc_gridsize_id, nc_grdcntrlat_id)
695    call netcdf_error_handler(ncstat)
696    ncstat = nf90_put_att(nc_grid_id, nc_grdcntrlat_id, 'units', 'degrees')
697    call netcdf_error_handler(ncstat)
698 
699    ! -----------------------------------------------------------------------------
700    ! - define grid center longitude array
701 
702    ncstat = nf90_def_var(nc_grid_id, 'grid_center_lon', NF90_DOUBLE, &
703                          nc_gridsize_id, nc_grdcntrlon_id)
704    call netcdf_error_handler(ncstat)
705    ncstat = nf90_put_att(nc_grid_id, nc_grdcntrlon_id, 'units', 'degrees')
706    call netcdf_error_handler(ncstat)
707 
708    ! -----------------------------------------------------------------------------
709    ! - define grid corner latitude array
710 
711    grid_dim_ids = (/ nc_gridcorn_id, nc_gridsize_id /)
712    ncstat = nf90_def_var(nc_grid_id, 'grid_corner_lat', NF90_DOUBLE, &
713                          grid_dim_ids, nc_grdcrnrlat_id)
714    call netcdf_error_handler(ncstat)
715    ncstat = nf90_put_att(nc_grid_id, nc_grdcrnrlat_id, 'units', 'degrees')
716    call netcdf_error_handler(ncstat)
717 
718    ! -----------------------------------------------------------------------------
719    ! - define grid corner longitude array
720 
721    ncstat = nf90_def_var(nc_grid_id, 'grid_corner_lon', NF90_DOUBLE, &
722                          grid_dim_ids, nc_grdcrnrlon_id)
723    call netcdf_error_handler(ncstat)
724    ncstat = nf90_put_att(nc_grid_id, nc_grdcrnrlon_id, 'units', 'degrees')
725    call netcdf_error_handler(ncstat)
726 
727    ! -----------------------------------------------------------------------------
728    ! end definition stage
729 
730    ncstat = nf90_enddef(nc_grid_id)
731    call netcdf_error_handler(ncstat)
732 
733    ! -----------------------------------------------------------------------------
734    ! write grid data
735 
736    ncstat = nf90_put_var(nc_grid_id, nc_griddims_id, grid_dims)
737    call netcdf_error_handler(ncstat)
738    ncstat = nf90_put_var(nc_grid_id, nc_grdimask_id, grid_imask)
739    call netcdf_error_handler(ncstat)
740    ncstat = nf90_put_var(nc_grid_id, nc_grdcntrlat_id, grid_center_lat)
741    call netcdf_error_handler(ncstat)
742    ncstat = nf90_put_var(nc_grid_id, nc_grdcntrlon_id, grid_center_lon)
743    call netcdf_error_handler(ncstat)
744    ncstat = nf90_put_var(nc_grid_id, nc_grdcrnrlat_id, grid_corner_lat)
745    call netcdf_error_handler(ncstat)
746    ncstat = nf90_put_var(nc_grid_id, nc_grdcrnrlon_id, grid_corner_lon)
747    call netcdf_error_handler(ncstat)
748 
749    ncstat = nf90_close(nc_grid_id)
750    call netcdf_error_handler(ncstat)
751 
752    DEALLOCATE( grid_imask, grid_center_lon, grid_center_lat, &
753                grid_corner_lon, grid_corner_lat )
754 
755 
756  END SUBROUTINE createSCRIPgrid
757 
758END MODULE scripgrid_mod
759
Note: See TracBrowser for help on using the repository browser.