1 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
2 | ! |
---|
3 | ! This module creates grid description files for input to the SCRIP code |
---|
4 | ! |
---|
5 | !----------------------------------------------------------------------- |
---|
6 | |
---|
7 | MODULE 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 | |
---|
47 | CONTAINS |
---|
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 | |
---|
758 | END MODULE scripgrid_mod |
---|
759 | |
---|