weight_gen  1.12.0
 All Files Functions
scrip.F90
Go to the documentation of this file.
1 
5 
19  program scrip
20 
21  implicit none
22 
23  character(len=128) :: outfile
24  character(len=20) :: title
25  character(len=5) :: idim_ch, jdim_ch, jdimp_ch
26  character(len=6) :: cres
27 
28  integer :: header_buffer_val = 16384
29  integer :: fsize=65536, inital = 0
30  integer :: error, ncid
31  integer :: i, j, idim, jdim, ijdim
32  integer :: jdimp, num_ch
33  integer :: dim_size, dim_corners, dim_rank
34  integer :: id_dims, id_center_lat, id_center_lon
35  integer :: id_imask, id_corner_lat, id_corner_lon
36  integer :: num_corners = 4
37  integer :: rank = 2
38  integer(kind=4), allocatable :: mask(:)
39 
40  real(kind=8) :: corner_lon_src
41  real(kind=8) :: dx_src, lat_edge
42  real(kind=8), allocatable :: lats(:,:), lons(:,:), dum1d(:)
43  real(kind=8), allocatable :: dum2d(:,:), latsp(:,:), lonsp(:,:)
44  real(kind=8), allocatable :: lats_corner(:,:,:), lons_corner(:,:,:)
45  real(kind=8), allocatable :: latsp_corner(:,:,:), lonsp_corner(:,:,:)
46  real(kind=8), allocatable :: slat(:), wlat(:)
47 
48  include "netcdf.inc"
49 
50  call getarg(1, cres)
51 
52  select case (trim(cres))
53  case ("c48","C48")
54  idim = 192 ! cres * 4
55  jdim = 94 ! cres * 2
56  jdimp = 96 ! include two rows for the poles
57  idim_ch = "192"
58  jdim_ch = "94"
59  jdimp_ch = "96"
60  case ("c96","C96")
61  idim = 384 ! cres * 4
62  jdim = 192 ! cres * 2
63  jdimp = 194 ! include two rows for the poles
64  idim_ch = "384"
65  jdim_ch = "192"
66  jdimp_ch = "194"
67  case ("c128","C128")
68  idim = 512 ! cres * 4
69  jdim = 256 ! cres * 2
70  jdimp = 258 ! include two rows for the poles
71  idim_ch = "512"
72  jdim_ch = "256"
73  jdimp_ch = "258"
74  case ("c192","C192")
75  idim = 768 ! cres * 4
76  jdim = 384 ! cres * 2
77  jdimp = 386 ! include two rows for the poles
78  idim_ch = "768"
79  jdim_ch = "384"
80  jdimp_ch = "386"
81  case ("c384","C384")
82  idim = 1536 ! cres * 4
83  jdim = 768 ! cres * 2
84  jdimp = 770 ! include two rows for the poles
85  idim_ch = "1536"
86  jdim_ch = "768"
87  jdimp_ch = "770"
88  case ("c768","C768")
89  idim = 3072 ! cres * 4
90  jdim = 1536 ! cres * 2
91  jdimp = 1538 ! include two rows for the poles
92  idim_ch = "3072"
93  jdim_ch = "1536"
94  jdimp_ch = "1538"
95  case ("c1152","C1152")
96  idim = 4608 ! cres * 4
97  jdim = 2304 ! cres * 2
98  jdimp = 2306 ! include two rows for the poles
99  idim_ch = "4608"
100  jdim_ch = "2304"
101  jdimp_ch = "2306"
102  case ("c3072","C3072")
103  idim = 12288 ! cres * 4
104  jdim = 6144 ! cres * 2
105  jdimp = 6146 ! include two rows for the poles
106  idim_ch = "12288"
107  jdim_ch = "6144"
108  jdimp_ch = "6146"
109  case default
110  print*,'- Resolution not supported ', trim(cres)
111  stop 3
112  end select
113 
114  corner_lon_src = 0.0
115  dx_src = 360.0 / float(idim)
116  ijdim = idim*jdim
117 
118  allocate(slat(jdim))
119  allocate(wlat(jdim))
120 
121  call splat(4, jdim, slat, wlat)
122 
123  allocate(lats(idim,jdim))
124  allocate(lats_corner(num_corners,idim,jdim))
125  allocate(lons(idim,jdim))
126  allocate(lons_corner(num_corners,idim,jdim))
127 
128  do j = 1, jdim
129  lats(:,j) = 90.0 - (acos(slat(j))* 180.0 / (4.*atan(1.)))
130  enddo
131 
132  deallocate(slat, wlat)
133 
134 !----------------------------------------------------------------
135 ! First, output file without poles.
136 !----------------------------------------------------------------
137 
138 !----------------------------------------------------------------
139 ! Set corners in counter-clockwise order
140 !
141 ! 2 1
142 !
143 ! C
144 !
145 ! 3 4
146 !----------------------------------------------------------------
147 
148  lats_corner(1,:,1) = 90.0
149  lats_corner(2,:,1) = 90.0
150 
151  lats_corner(3,:,jdim) = -90.0
152  lats_corner(4,:,jdim) = -90.0
153 
154  do j = 1, jdim - 1
155  lat_edge = (lats(1,j) + lats(1,j+1)) / 2.0
156  lats_corner(3,:,j) = lat_edge
157  lats_corner(4,:,j) = lat_edge
158  lats_corner(1,:,j+1) = lat_edge
159  lats_corner(2,:,j+1) = lat_edge
160  enddo
161 
162  do i = 1, idim
163  lons(i,:) = corner_lon_src + float(i-1)*dx_src
164  lons_corner(1,i,:) = lons(i,:) + (dx_src*0.5)
165  lons_corner(2,i,:) = lons(i,:) - (dx_src*0.5)
166  lons_corner(3,i,:) = lons(i,:) - (dx_src*0.5)
167  lons_corner(4,i,:) = lons(i,:) + (dx_src*0.5)
168  enddo
169 
170  i = 1
171  j = 1
172  print*,'center ',lats(i,j),lons(i,j)
173  print*,'corner 1 ',lats_corner(1,i,j),lons_corner(1,i,j)
174  print*,'corner 2 ',lats_corner(2,i,j),lons_corner(2,i,j)
175  print*,'corner 3 ',lats_corner(3,i,j),lons_corner(3,i,j)
176  print*,'corner 4 ',lats_corner(4,i,j),lons_corner(4,i,j)
177 
178  i = 1
179  j = 2
180  print*,'center ',lats(i,j),lons(i,j)
181  print*,'corner 1 ',lats_corner(1,i,j),lons_corner(1,i,j)
182  print*,'corner 2 ',lats_corner(2,i,j),lons_corner(2,i,j)
183  print*,'corner 3 ',lats_corner(3,i,j),lons_corner(3,i,j)
184  print*,'corner 4 ',lats_corner(4,i,j),lons_corner(4,i,j)
185 
186  i = 1
187  j = jdim
188  print*,'center ',lats(i,j),lons(i,j)
189  print*,'corner 1 ',lats_corner(1,i,j),lons_corner(1,i,j)
190  print*,'corner 2 ',lats_corner(2,i,j),lons_corner(2,i,j)
191  print*,'corner 3 ',lats_corner(3,i,j),lons_corner(3,i,j)
192  print*,'corner 4 ',lats_corner(4,i,j),lons_corner(4,i,j)
193 
194  i = 1
195  j = jdim-1
196  print*,'center ',lats(i,j),lons(i,j)
197  print*,'corner 1 ',lats_corner(1,i,j),lons_corner(1,i,j)
198  print*,'corner 2 ',lats_corner(2,i,j),lons_corner(2,i,j)
199  print*,'corner 3 ',lats_corner(3,i,j),lons_corner(3,i,j)
200  print*,'corner 4 ',lats_corner(4,i,j),lons_corner(4,i,j)
201 
202  allocate(mask(ijdim))
203  mask = 1
204 
205 ! output file without pole.
206 
207  outfile = " "
208  outfile = "./gaussian." // trim(idim_ch) // "." // trim(jdim_ch) // ".nc"
209  title = " "
210  title = "gaussian." // trim(idim_ch) // "." // trim(jdim_ch)
211 
212 !--- open the file
213  error = nf__create(outfile, ior(nf_netcdf4,nf_classic_model), inital, fsize, ncid)
214  print*, 'error after open ', error
215 
216 !--- define dimension
217  error = nf_def_dim(ncid, 'grid_size', ijdim, dim_size)
218  error = nf_def_dim(ncid, 'grid_corners', num_corners, dim_corners)
219  error = nf_def_dim(ncid, 'grid_rank', rank, dim_rank)
220 
221 !--- define field
222  error = nf_def_var(ncid, 'grid_dims', nf_int, 1, (/dim_rank/), id_dims)
223  error = nf_def_var(ncid, 'grid_center_lat', nf_double, 1, (/dim_size/), id_center_lat)
224  error = nf_put_att_text(ncid, id_center_lat, "units", 7, "degrees")
225  error = nf_def_var(ncid, 'grid_center_lon', nf_double, 1, (/dim_size/), id_center_lon)
226  error = nf_put_att_text(ncid, id_center_lon, "units", 7, "degrees")
227  error = nf_def_var(ncid, 'grid_imask', nf_int, 1, (/dim_size/), id_imask)
228  error = nf_put_att_text(ncid, id_imask, "units", 8, "unitless")
229  error = nf_def_var(ncid, 'grid_corner_lat', nf_double, 2, (/dim_corners,dim_size/), id_corner_lat)
230  error = nf_put_att_text(ncid, id_corner_lat, "units", 7, "degrees")
231  error = nf_def_var(ncid, 'grid_corner_lon', nf_double, 2, (/dim_corners,dim_size/), id_corner_lon)
232  error = nf_put_att_text(ncid, id_corner_lon, "units", 7, "degrees")
233  num_ch = len_trim(title)
234  error = nf_put_att_text(ncid, nf_global, "title", num_ch, trim(title))
235  error = nf__enddef(ncid, header_buffer_val,4,0,4)
236 
237 !--- set fields
238  error = nf_put_var_int( ncid, id_dims, (/idim,jdim/))
239 
240  allocate(dum1d(ijdim))
241  dum1d = reshape(lats, (/ijdim/))
242  error = nf_put_var_double( ncid, id_center_lat, dum1d)
243  dum1d = reshape(lons, (/ijdim/))
244  error = nf_put_var_double( ncid, id_center_lon, dum1d)
245  deallocate(dum1d)
246 
247  error = nf_put_var_int( ncid, id_imask, mask)
248  deallocate(mask)
249 
250  allocate(dum2d(num_corners,ijdim))
251  dum2d = reshape(lats_corner, (/num_corners,ijdim/))
252  error = nf_put_var_double( ncid, id_corner_lat, dum2d)
253 
254  dum2d = reshape(lons_corner, (/num_corners,ijdim/))
255  error = nf_put_var_double( ncid, id_corner_lon, dum2d)
256  deallocate(dum2d)
257 
258  error = nf_close(ncid)
259 
260 !----------------------------------------------------------------
261 ! output file with poles.
262 !----------------------------------------------------------------
263 
264  outfile = " "
265  outfile = "./gaussian." // trim(idim_ch) // "." // trim(jdimp_ch) // ".nc"
266  title = " "
267  title = "gaussian." // trim(idim_ch) // "." // trim(jdimp_ch)
268 
269  ijdim = idim*jdimp
270 
271  allocate(latsp(idim,jdimp))
272  allocate(lonsp(idim,jdimp))
273 
274  do j = 2, jdim+1
275  latsp(:,j) = lats(:,j-1)
276  lonsp(:,j) = lons(:,j-1)
277  enddo
278 
279  latsp(:,1) = 90.0_8
280  lonsp(:,1) = 0.0_8
281 
282  latsp(:,jdimp) = -90.0_8
283  lonsp(:,jdimp) = 0.0_8
284 
285  deallocate(lats, lons)
286 
287  allocate(latsp_corner(num_corners,idim,jdimp))
288  allocate(lonsp_corner(num_corners,idim,jdimp))
289 
290  latsp_corner(:,:,1) = 89.5_8
291  latsp_corner(:,:,jdimp) = -89.5_8
292 
293  lonsp_corner(1,:,1) = 0.0_8
294  lonsp_corner(2,:,1) = 90.0_8
295  lonsp_corner(3,:,1) = 180.0_8
296  lonsp_corner(4,:,1) = 270.0_8
297 
298  lonsp_corner(1,:,jdimp) = 0.0_8
299  lonsp_corner(2,:,jdimp) = 90.0_8
300  lonsp_corner(3,:,jdimp) = 180.0_8
301  lonsp_corner(4,:,jdimp) = 270.0_8
302 
303  do j = 2, jdim+1
304  latsp_corner(:,:,j) = lats_corner(:,:,j-1)
305  lonsp_corner(:,:,j) = lons_corner(:,:,j-1)
306  enddo
307 
308  deallocate(lats_corner, lons_corner)
309 
310 !--- open the file
311  error = nf__create(outfile, ior(nf_netcdf4,nf_classic_model), inital, fsize, ncid)
312  print*, 'error after open ', error
313 
314 !--- define dimension
315  error = nf_def_dim(ncid, 'grid_size', ijdim, dim_size)
316  error = nf_def_dim(ncid, 'grid_corners', num_corners, dim_corners)
317  error = nf_def_dim(ncid, 'grid_rank', rank, dim_rank)
318 
319 !--- define field
320  error = nf_def_var(ncid, 'grid_dims', nf_int, 1, (/dim_rank/), id_dims)
321  error = nf_def_var(ncid, 'grid_center_lat', nf_double, 1, (/dim_size/), id_center_lat)
322  error = nf_put_att_text(ncid, id_center_lat, "units", 7, "degrees")
323  error = nf_def_var(ncid, 'grid_center_lon', nf_double, 1, (/dim_size/), id_center_lon)
324  error = nf_put_att_text(ncid, id_center_lon, "units", 7, "degrees")
325  error = nf_def_var(ncid, 'grid_imask', nf_int, 1, (/dim_size/), id_imask)
326  error = nf_put_att_text(ncid, id_imask, "units", 8, "unitless")
327  error = nf_def_var(ncid, 'grid_corner_lat', nf_double, 2, (/dim_corners,dim_size/), id_corner_lat)
328  error = nf_put_att_text(ncid, id_corner_lat, "units", 7, "degrees")
329  error = nf_def_var(ncid, 'grid_corner_lon', nf_double, 2, (/dim_corners,dim_size/), id_corner_lon)
330  error = nf_put_att_text(ncid, id_corner_lon, "units", 7, "degrees")
331  num_ch = len_trim(title)
332  error = nf_put_att_text(ncid, nf_global, "title", num_ch, trim(title))
333  error = nf__enddef(ncid, header_buffer_val,4,0,4)
334 
335 !--- set fields
336  error = nf_put_var_int( ncid, id_dims, (/idim,jdimp/))
337 
338  allocate(dum1d(ijdim))
339  dum1d = reshape(latsp, (/ijdim/))
340  error = nf_put_var_double( ncid, id_center_lat, dum1d)
341  dum1d = reshape(lonsp, (/ijdim/))
342  error = nf_put_var_double( ncid, id_center_lon, dum1d)
343  deallocate(dum1d)
344 
345  allocate(mask(ijdim))
346  mask = 1
347  error = nf_put_var_int( ncid, id_imask, mask)
348  deallocate(mask)
349 
350  allocate(dum2d(num_corners,ijdim))
351  dum2d = reshape(latsp_corner, (/num_corners,ijdim/))
352  print*,'lat corner check ',maxval(dum2d),minval(dum2d)
353  error = nf_put_var_double( ncid, id_corner_lat, dum2d)
354  deallocate(latsp_corner)
355 
356  dum2d = reshape(lonsp_corner, (/num_corners,ijdim/))
357  error = nf_put_var_double( ncid, id_corner_lon, dum2d)
358  deallocate(dum2d, lonsp_corner)
359 
360  error = nf_close(ncid)
361 
362  print*,'- DONE.'
363 
364  end program scrip
scrip
program scrip
Create "scrip" files that describes a gaussian grid.
Definition: scrip.F90:19