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.
scripgrid_mod.F90 in branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/src – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/src/scripgrid_mod.F90 @ 2352

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