weight_gen 1.14.0
Loading...
Searching...
No Matches
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