source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/SCRIP1.4/grids/convertPOPT.f @ 2352

Last change on this file since 2352 was 2352, checked in by sga, 11 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: 14.8 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     This file converts a POP grid.dat file to a remapping grid file
4!     in netCDF format.
5!
6!-----------------------------------------------------------------------
7!
8!     CVS:$Id: convertPOPT.f,v 1.4 2001/08/21 21:22:56 pwjones Exp $
9!
10!     Copyright (c) 1997, 1998 the Regents of the University of
11!       California.
12!
13!     Unless otherwise indicated, this software has been authored
14!     by an employee or employees of the University of California,
15!     operator of the Los Alamos National Laboratory under Contract
16!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
17!     Government has rights to use, reproduce, and distribute this
18!     software.  The public may copy and use this software without
19!     charge, provided that this Notice and any statement of authorship
20!     are reproduced on all copies.  Neither the Government nor the
21!     University makes any warranty, express or implied, or assumes
22!     any liability or responsibility for the use of this software.
23!
24!***********************************************************************
25
26      program convertPOPT
27
28!-----------------------------------------------------------------------
29!
30!     This file converts a POP grid.dat file to a remapping grid file.
31!
32!-----------------------------------------------------------------------
33
34      use kinds_mod
35      use constants
36      use iounits
37      use netcdf_mod
38
39      implicit none
40
41!-----------------------------------------------------------------------
42!
43!     variables that describe the grid
44!       4/3       nx = 192, ny = 128
45!       2/3 (mod) nx = 384, ny = 288
46!       x3p Greenland DP nx = 100, ny = 116
47!       x2p Greenland DP nx = 160, ny = 192
48!       x1p Greenland DP nx = 320, ny = 384
49!
50!-----------------------------------------------------------------------
51
52      integer (kind=int_kind), parameter ::
53     &             nx = 320, ny = 384,
54     &             grid_size = nx*ny,
55     &             grid_rank = 2,
56     &             grid_corners = 4
57
58      integer (kind=int_kind), dimension(2) ::
59     &             grid_dims   ! size of each dimension
60
61      character(char_len), parameter :: 
62     &    grid_name = 'Greenland DP x1p',
63     &    grid_file_in  = '/scratch/pwjones/grid.320x384.da',
64     &    grid_topo_in  = '/scratch/pwjones/kmt.320x384.da',
65     &    grid_file_out = '/scratch/pwjones/Greenland_DP_x1p.nc'
66
67      real (kind=dbl_kind), parameter ::
68     &  radius    = 6370.0e5_dbl_kind       ! radius of Earth (cm)
69     &, area_norm = one/(radius*radius)
70
71!-----------------------------------------------------------------------
72!
73!     grid coordinates and masks
74!
75!-----------------------------------------------------------------------
76
77      integer (kind=int_kind), dimension(grid_size) ::
78     &             grid_imask
79
80      real (kind=dbl_kind), dimension(grid_size) ::
81     &             grid_area      ,  ! area as computed in POP
82     &             grid_center_lat,  ! lat/lon coordinates for
83     &             grid_center_lon   ! each grid center in radians
84
85      real (kind=dbl_kind), dimension(grid_corners,grid_size) ::
86     &             grid_corner_lat,  ! lat/lon coordinates for
87     &             grid_corner_lon   ! each grid corner in radians
88
89      real (kind=dbl_kind), dimension(nx,ny) ::
90     &             HTN, HTE          ! T-cell grid lengths
91
92!-----------------------------------------------------------------------
93!
94!     other local variables
95!
96!-----------------------------------------------------------------------
97
98      integer (kind=int_kind) :: i, j, n, iunit, ocn_add, im1, jm1
99
100      integer (kind=int_kind) ::
101     &        ncstat,            ! general netCDF status variable
102     &        nc_grid_id,        ! netCDF grid dataset id
103     &        nc_gridsize_id,    ! netCDF grid size dim id
104     &        nc_gridcorn_id,    ! netCDF grid corner dim id
105     &        nc_gridrank_id,    ! netCDF grid rank dim id
106     &        nc_griddims_id,    ! netCDF grid dimensions id
107     &        nc_grdcntrlat_id,  ! netCDF grid center lat id
108     &        nc_grdcntrlon_id,  ! netCDF grid center lon id
109     &        nc_grdimask_id,    ! netCDF grid mask id
110     &        nc_gridarea_id,    ! netCDF grid area id
111     &        nc_grdcrnrlat_id,  ! netCDF grid corner lat id
112     &        nc_grdcrnrlon_id   ! netCDF grid corner lon id
113
114      integer (kind=int_kind), dimension(2) ::
115     &        nc_dims2_id        ! netCDF dim id array for 2-d arrays
116
117      real (kind=dbl_kind) :: tmplon, dxt, dyt
118
119!-----------------------------------------------------------------------
120!
121!     read in grid info
122!     lat/lon info is on velocity points which correspond
123!     to the NE corner (in logical space) of the grid cell.
124!
125!-----------------------------------------------------------------------
126
127      call get_unit(iunit)
128      open(unit=iunit, file=grid_topo_in, status='old', 
129     &     form='unformatted', access='direct', recl=grid_size*4)
130      read (unit=iunit,rec=1) grid_imask
131      call release_unit(iunit)
132
133      call get_unit(iunit)
134      open(unit=iunit, file=grid_file_in, status='old', 
135     &     form='unformatted', access='direct', recl=grid_size*8)
136      read (unit=iunit, rec=1) grid_corner_lat(3,:)
137      read (unit=iunit, rec=2) grid_corner_lon(3,:)
138      read (unit=iunit, rec=3) HTN
139      read (unit=iunit, rec=4) HTE
140      call release_unit(iunit)
141
142      grid_dims(1) = nx
143      grid_dims(2) = ny
144
145!-----------------------------------------------------------------------
146!
147!     convert KMT field to integer grid mask
148!
149!-----------------------------------------------------------------------
150
151      grid_imask = min(grid_imask, 1)
152
153!-----------------------------------------------------------------------
154!
155!     compute remaining corners
156!
157!-----------------------------------------------------------------------
158
159      do j=1,ny
160        do i=1,nx
161          ocn_add = (j-1)*nx + i
162          if (i .ne. 1) then
163            im1 = ocn_add - 1
164          else
165            im1 = ocn_add + nx - 1
166          endif
167
168          grid_corner_lat(4,ocn_add) = grid_corner_lat(3,im1)
169          grid_corner_lon(4,ocn_add) = grid_corner_lon(3,im1)
170        end do
171      end do
172
173      do j=2,ny
174        do i=1,nx
175          ocn_add = (j-1)*nx + i
176          jm1 = (j-2)*nx + i
177
178          grid_corner_lat(2,ocn_add) = grid_corner_lat(3,jm1)
179          grid_corner_lat(1,ocn_add) = grid_corner_lat(4,jm1)
180
181          grid_corner_lon(2,ocn_add) = grid_corner_lon(3,jm1)
182          grid_corner_lon(1,ocn_add) = grid_corner_lon(4,jm1)
183        end do
184      end do
185
186!-----------------------------------------------------------------------
187!
188!     mock up the lower row boundaries
189!
190!-----------------------------------------------------------------------
191
192      do i=1,nx
193        grid_corner_lat(1,i) = -pih + tiny
194        grid_corner_lat(2,i) = -pih + tiny
195
196        grid_corner_lon(1,i) = grid_corner_lon(4,i)
197        grid_corner_lon(2,i) = grid_corner_lon(3,i)
198      end do
199
200!-----------------------------------------------------------------------
201!
202!     correct for 0,2pi longitude crossings
203!
204!-----------------------------------------------------------------------
205
206      do ocn_add=1,grid_size
207        if (grid_corner_lon(1,ocn_add) > pi2) 
208     &      grid_corner_lon(1,ocn_add) = 
209     &      grid_corner_lon(1,ocn_add) - pi2
210        if (grid_corner_lon(1,ocn_add) < 0.0) 
211     &      grid_corner_lon(1,ocn_add) = 
212     &      grid_corner_lon(1,ocn_add) + pi2
213        do n=2,grid_corners
214          tmplon = grid_corner_lon(n  ,ocn_add) - 
215     &             grid_corner_lon(n-1,ocn_add) 
216          if (tmplon < -three*pih) grid_corner_lon(n,ocn_add) = 
217     &                             grid_corner_lon(n,ocn_add) + pi2
218          if (tmplon >  three*pih) grid_corner_lon(n,ocn_add) = 
219     &                             grid_corner_lon(n,ocn_add) - pi2
220        end do
221      end do
222
223!-----------------------------------------------------------------------
224!
225!     compute ocean cell centers by averaging corner values
226!
227!-----------------------------------------------------------------------
228
229      do ocn_add=1,grid_size
230        grid_center_lat(ocn_add) = grid_corner_lat(1,ocn_add)
231        grid_center_lon(ocn_add) = grid_corner_lon(1,ocn_add)
232        do n=2,grid_corners
233          grid_center_lat(ocn_add) = grid_center_lat(ocn_add) + 
234     &                               grid_corner_lat(n,ocn_add)
235          grid_center_lon(ocn_add) = grid_center_lon(ocn_add) + 
236     &                               grid_corner_lon(n,ocn_add)
237        end do
238        grid_center_lat(ocn_add) = grid_center_lat(ocn_add)/
239     &                                  float(grid_corners)
240        grid_center_lon(ocn_add) = grid_center_lon(ocn_add)/
241     &                                  float(grid_corners)
242        if (grid_center_lon(ocn_add) > pi2) 
243     &      grid_center_lon(ocn_add) = grid_center_lon(ocn_add) - pi2
244        if (grid_center_lon(ocn_add) < 0.0) 
245     &      grid_center_lon(ocn_add) = grid_center_lon(ocn_add) + pi2
246      end do
247
248!-----------------------------------------------------------------------
249!
250!     compute cell areas in same way as POP
251!
252!-----------------------------------------------------------------------
253
254      n = 0
255      do j=1,ny
256        if (j > 1) then
257          jm1 = j-1
258        else
259          jm1 = 1
260        endif
261        do i=1,nx
262          if (i > 1) then
263            im1 = i-1
264          else
265            im1 = nx
266          endif
267
268          n = n+1
269
270          dxt = half*(HTN(i,j) + HTN(i,jm1))
271          dyt = half*(HTE(i,j) + HTE(im1,j))
272          if (dxt == zero) dxt=one
273          if (dyt == zero) dyt=one
274
275          grid_area(n) = dxt*dyt*area_norm
276        end do
277      end do
278
279!-----------------------------------------------------------------------
280!
281!     set up attributes for netCDF file
282!
283!-----------------------------------------------------------------------
284
285      !***
286      !*** create netCDF dataset for this grid
287      !***
288
289      ncstat = nf_create (grid_file_out, NF_CLOBBER,
290     &                    nc_grid_id)
291      call netcdf_error_handler(ncstat)
292
293      ncstat = nf_put_att_text (nc_grid_id, NF_GLOBAL, 'title',
294     &                          len_trim(grid_name), grid_name)
295      call netcdf_error_handler(ncstat)
296
297      !***
298      !*** define grid size dimension
299      !***
300
301      ncstat = nf_def_dim (nc_grid_id, 'grid_size', grid_size, 
302     &                     nc_gridsize_id)
303      call netcdf_error_handler(ncstat)
304
305      !***
306      !*** define grid rank dimension
307      !***
308
309      ncstat = nf_def_dim (nc_grid_id, 'grid_rank', grid_rank, 
310     &                     nc_gridrank_id)
311      call netcdf_error_handler(ncstat)
312
313      !***
314      !*** define grid corner dimension
315      !***
316
317      ncstat = nf_def_dim (nc_grid_id, 'grid_corners', grid_corners, 
318     &                     nc_gridcorn_id)
319      call netcdf_error_handler(ncstat)
320
321      !***
322      !*** define grid dim size array
323      !***
324
325      ncstat = nf_def_var (nc_grid_id, 'grid_dims', NF_INT,
326     &                     1, nc_gridrank_id, nc_griddims_id)
327      call netcdf_error_handler(ncstat)
328
329      !***
330      !*** define grid center latitude array
331      !***
332
333      ncstat = nf_def_var (nc_grid_id, 'grid_center_lat', NF_DOUBLE,
334     &                     1, nc_gridsize_id, nc_grdcntrlat_id)
335      call netcdf_error_handler(ncstat)
336
337      ncstat = nf_put_att_text (nc_grid_id, nc_grdcntrlat_id, 'units',
338     &                          7, 'radians')
339      call netcdf_error_handler(ncstat)
340
341      !***
342      !*** define grid center longitude array
343      !***
344
345      ncstat = nf_def_var (nc_grid_id, 'grid_center_lon', NF_DOUBLE,
346     &                     1, nc_gridsize_id, nc_grdcntrlon_id)
347      call netcdf_error_handler(ncstat)
348
349      ncstat = nf_put_att_text (nc_grid_id, nc_grdcntrlon_id, 'units',
350     &                          7, 'radians')
351      call netcdf_error_handler(ncstat)
352
353      !***
354      !*** define grid area array
355      !***
356
357      ncstat = nf_def_var (nc_grid_id, 'grid_area', NF_DOUBLE,
358     &                     1, nc_gridsize_id, nc_gridarea_id)
359      call netcdf_error_handler(ncstat)
360
361      !***
362      !*** define grid mask
363      !***
364
365      ncstat = nf_def_var (nc_grid_id, 'grid_imask', NF_INT,
366     &                     1, nc_gridsize_id, nc_grdimask_id)
367      call netcdf_error_handler(ncstat)
368
369      ncstat = nf_put_att_text (nc_grid_id, nc_grdimask_id, 'units',
370     &                          8, 'unitless')
371      call netcdf_error_handler(ncstat)
372
373      !***
374      !*** define grid corner latitude array
375      !***
376
377      nc_dims2_id(1) = nc_gridcorn_id
378      nc_dims2_id(2) = nc_gridsize_id
379
380      ncstat = nf_def_var (nc_grid_id, 'grid_corner_lat', NF_DOUBLE,
381     &                     2, nc_dims2_id, nc_grdcrnrlat_id)
382      call netcdf_error_handler(ncstat)
383
384      ncstat = nf_put_att_text (nc_grid_id, nc_grdcrnrlat_id, 'units',
385     &                          7, 'radians')
386      call netcdf_error_handler(ncstat)
387
388      !***
389      !*** define grid corner longitude array
390      !***
391
392      ncstat = nf_def_var (nc_grid_id, 'grid_corner_lon', NF_DOUBLE,
393     &                     2, nc_dims2_id, nc_grdcrnrlon_id)
394      call netcdf_error_handler(ncstat)
395
396      ncstat = nf_put_att_text (nc_grid_id, nc_grdcrnrlon_id, 'units',
397     &                          7, 'radians')
398      call netcdf_error_handler(ncstat)
399
400      !***
401      !*** end definition stage
402      !***
403
404      ncstat = nf_enddef(nc_grid_id)
405      call netcdf_error_handler(ncstat)
406
407!-----------------------------------------------------------------------
408!
409!     write grid data
410!
411!-----------------------------------------------------------------------
412
413      ncstat = nf_put_var_int(nc_grid_id, nc_griddims_id, grid_dims)
414      call netcdf_error_handler(ncstat)
415
416      ncstat = nf_put_var_double(nc_grid_id, nc_grdcntrlat_id, 
417     &                           grid_center_lat)
418      ncstat = nf_put_var_int(nc_grid_id, nc_grdimask_id, grid_imask)
419      call netcdf_error_handler(ncstat)
420
421      ncstat = nf_put_var_double(nc_grid_id, nc_gridarea_id, 
422     &                           grid_area)
423      call netcdf_error_handler(ncstat)
424
425      ncstat = nf_put_var_double(nc_grid_id, nc_grdcntrlat_id, 
426     &                           grid_center_lat)
427      call netcdf_error_handler(ncstat)
428
429      ncstat = nf_put_var_double(nc_grid_id, nc_grdcntrlon_id, 
430     &                           grid_center_lon)
431      call netcdf_error_handler(ncstat)
432
433      ncstat = nf_put_var_double(nc_grid_id, nc_grdcrnrlat_id, 
434     &                           grid_corner_lat)
435      call netcdf_error_handler(ncstat)
436
437      ncstat = nf_put_var_double(nc_grid_id, nc_grdcrnrlon_id, 
438     &                           grid_corner_lon)
439      call netcdf_error_handler(ncstat)
440
441      ncstat = nf_close(nc_grid_id)
442
443!***********************************************************************
444
445      end program convertPOPT
446
447!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.