source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/SCRIP1.4/grids/convertgauss.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: 15.2 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     This program creates a remapping grid file for Gaussian lat/lon
4!     grids (for spectral transform codes).
5!
6!-----------------------------------------------------------------------
7!
8!     CVS:$Id: convertgauss.f,v 1.3 2000/04/19 22:05:57 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 convert_gauss
27
28!-----------------------------------------------------------------------
29!
30!     This file creates a remapping grid file for a Gaussian grid
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!
45!     T42: nx=128 ny=64
46!     T62: nx=192 ny=94
47!
48!-----------------------------------------------------------------------
49
50      integer (kind=int_kind), parameter ::
51     &             nx = 192, ny = 94,
52     &             grid_size = nx*ny,
53     &             grid_rank = 2,
54     &             grid_corners = 4
55
56      character(char_len), parameter :: 
57     &             grid_name = 'T62 Gaussian Grid',
58     &             grid_file_out = 'remap_grid_T62.nc'
59
60      integer (kind=int_kind), dimension(grid_rank) ::
61     &             grid_dims
62
63!-----------------------------------------------------------------------
64!
65!     grid coordinates and masks
66!
67!-----------------------------------------------------------------------
68
69      integer (kind=int_kind), dimension(grid_size) ::
70     &             grid_imask
71
72      real (kind=dbl_kind), dimension(grid_size) ::
73     &             grid_center_lat,  ! lat/lon coordinates for
74     &             grid_center_lon   ! each grid center in degrees
75
76      real (kind=dbl_kind), dimension(grid_corners,grid_size) ::
77     &             grid_corner_lat,  ! lat/lon coordinates for
78     &             grid_corner_lon   ! each grid corner in degrees
79
80!-----------------------------------------------------------------------
81!
82!     other local variables
83!
84!-----------------------------------------------------------------------
85
86      integer (kind=int_kind) :: i, j, iunit, atm_add
87
88      integer (kind=int_kind) ::
89     &        ncstat,            ! general netCDF status variable
90     &        nc_grid_id,        ! netCDF grid dataset id
91     &        nc_gridsize_id,    ! netCDF grid size dim id
92     &        nc_gridcorn_id,    ! netCDF grid corner dim id
93     &        nc_gridrank_id,    ! netCDF grid rank dim id
94     &        nc_griddims_id,    ! netCDF grid dimension size id
95     &        nc_grdcntrlat_id,  ! netCDF grid center lat id
96     &        nc_grdcntrlon_id,  ! netCDF grid center lon id
97     &        nc_grdimask_id,    ! netCDF grid mask id
98     &        nc_grdcrnrlat_id,  ! netCDF grid corner lat id
99     &        nc_grdcrnrlon_id   ! netCDF grid corner lon id
100
101      integer (kind=int_kind), dimension(2) ::
102     &        nc_dims2_id        ! netCDF dim id array for 2-d arrays
103
104      real (kind=dbl_kind) :: dlon, minlon, maxlon, centerlon,
105     &                              minlat, maxlat, centerlat
106
107      real (kind=dbl_kind), dimension(ny) :: gauss_root, gauss_wgt
108
109!-----------------------------------------------------------------------
110!
111!     compute longitudes of cell centers and corners.  set up alon
112!     array for search routine.
113!
114!-----------------------------------------------------------------------
115
116      grid_dims(1) = nx
117      grid_dims(2) = ny
118
119      dlon = 360./nx
120
121      do i=1,nx
122
123        centerlon = (i-1)*dlon
124        minlon = centerlon - half*dlon
125        maxlon = centerlon + half*dlon
126
127        do j=1,ny
128          atm_add = (j-1)*nx + i
129
130          grid_center_lon(atm_add  ) = centerlon
131          grid_corner_lon(1,atm_add) = minlon
132          grid_corner_lon(2,atm_add) = maxlon
133          grid_corner_lon(3,atm_add) = maxlon
134          grid_corner_lon(4,atm_add) = minlon
135        end do
136
137      end do
138
139!-----------------------------------------------------------------------
140!
141!     compute Gaussian latitudes and store in gauss_wgt.
142!
143!-----------------------------------------------------------------------
144
145      call gquad(ny, gauss_root, gauss_wgt)
146      do j=1,ny
147        gauss_wgt(j) = pih - gauss_root(ny+1-j)
148      end do
149
150!-----------------------------------------------------------------------
151!
152!     compute latitudes at cell centers and corners.  set up alat
153!     array for search routine.
154!
155!-----------------------------------------------------------------------
156
157      do j=1,ny
158        centerlat = gauss_wgt(j)
159
160        if (j .eq. 1) then
161          minlat = -pih
162        else
163          minlat = ATAN((COS(gauss_wgt(j-1)) - 
164     &                   COS(gauss_wgt(j  )))/
165     &                  (SIN(gauss_wgt(j  )) - 
166     &                   SIN(gauss_wgt(j-1))))
167        endif
168
169        if (j .eq. ny) then
170          maxlat = pih
171        else
172          maxlat = ATAN((COS(gauss_wgt(j  )) - 
173     &                   COS(gauss_wgt(j+1)))/
174     &                  (SIN(gauss_wgt(j+1)) - 
175     &                   SIN(gauss_wgt(j  ))))
176        endif
177
178        do i=1,nx
179          atm_add = (j-1)*nx + i
180          grid_center_lat(atm_add  ) = centerlat*360./pi2
181          grid_corner_lat(1,atm_add) = minlat*360./pi2
182          grid_corner_lat(2,atm_add) = minlat*360./pi2
183          grid_corner_lat(3,atm_add) = maxlat*360./pi2
184          grid_corner_lat(4,atm_add) = maxlat*360./pi2
185        end do
186
187      end do
188
189!-----------------------------------------------------------------------
190!
191!     define mask
192!
193!-----------------------------------------------------------------------
194
195      grid_imask = 1
196
197!-----------------------------------------------------------------------
198!
199!     set up attributes for netCDF file
200!
201!-----------------------------------------------------------------------
202
203      !***
204      !*** create netCDF dataset for this grid
205      !***
206
207      ncstat = nf_create (grid_file_out, NF_CLOBBER,
208     &                    nc_grid_id)
209      call netcdf_error_handler(ncstat)
210
211      ncstat = nf_put_att_text (nc_grid_id, NF_GLOBAL, 'title',
212     &                          len_trim(grid_name), grid_name)
213      call netcdf_error_handler(ncstat)
214
215      !***
216      !*** define grid size dimension
217      !***
218
219      ncstat = nf_def_dim (nc_grid_id, 'grid_size', grid_size, 
220     &                     nc_gridsize_id)
221      call netcdf_error_handler(ncstat)
222
223      !***
224      !*** define grid corner dimension
225      !***
226
227      ncstat = nf_def_dim (nc_grid_id, 'grid_corners', grid_corners, 
228     &                     nc_gridcorn_id)
229      call netcdf_error_handler(ncstat)
230
231      !***
232      !*** define grid rank dimension
233      !***
234
235      ncstat = nf_def_dim (nc_grid_id, 'grid_rank', grid_rank, 
236     &                     nc_gridrank_id)
237      call netcdf_error_handler(ncstat)
238
239      !***
240      !*** define grid dimension size array
241      !***
242
243      ncstat = nf_def_var (nc_grid_id, 'grid_dims', NF_INT,
244     &                     1, nc_gridrank_id, nc_griddims_id)
245      call netcdf_error_handler(ncstat)
246
247      !***
248      !*** define grid center latitude array
249      !***
250
251      ncstat = nf_def_var (nc_grid_id, 'grid_center_lat', NF_DOUBLE,
252     &                     1, nc_gridsize_id, nc_grdcntrlat_id)
253      call netcdf_error_handler(ncstat)
254
255      ncstat = nf_put_att_text (nc_grid_id, nc_grdcntrlat_id, 'units',
256     &                          7, 'degrees')
257      call netcdf_error_handler(ncstat)
258
259      !***
260      !*** define grid center longitude array
261      !***
262
263      ncstat = nf_def_var (nc_grid_id, 'grid_center_lon', NF_DOUBLE,
264     &                     1, nc_gridsize_id, nc_grdcntrlon_id)
265      call netcdf_error_handler(ncstat)
266
267      ncstat = nf_put_att_text (nc_grid_id, nc_grdcntrlon_id, 'units',
268     &                          7, 'degrees')
269      call netcdf_error_handler(ncstat)
270
271      !***
272      !*** define grid mask
273      !***
274
275      ncstat = nf_def_var (nc_grid_id, 'grid_imask', NF_INT,
276     &                     1, nc_gridsize_id, nc_grdimask_id)
277      call netcdf_error_handler(ncstat)
278
279      ncstat = nf_put_att_text (nc_grid_id, nc_grdimask_id, 'units',
280     &                          8, 'unitless')
281      call netcdf_error_handler(ncstat)
282
283      !***
284      !*** define grid corner latitude array
285      !***
286
287      nc_dims2_id(1) = nc_gridcorn_id
288      nc_dims2_id(2) = nc_gridsize_id
289
290      ncstat = nf_def_var (nc_grid_id, 'grid_corner_lat', NF_DOUBLE,
291     &                     2, nc_dims2_id, nc_grdcrnrlat_id)
292      call netcdf_error_handler(ncstat)
293
294      ncstat = nf_put_att_text (nc_grid_id, nc_grdcrnrlat_id, 'units',
295     &                          7, 'degrees')
296      call netcdf_error_handler(ncstat)
297
298      !***
299      !*** define grid corner longitude array
300      !***
301
302      ncstat = nf_def_var (nc_grid_id, 'grid_corner_lon', NF_DOUBLE,
303     &                     2, nc_dims2_id, nc_grdcrnrlon_id)
304      call netcdf_error_handler(ncstat)
305
306      ncstat = nf_put_att_text (nc_grid_id, nc_grdcrnrlon_id, 'units',
307     &                          7, 'degrees')
308      call netcdf_error_handler(ncstat)
309
310      !***
311      !*** end definition stage
312      !***
313
314      ncstat = nf_enddef(nc_grid_id)
315      call netcdf_error_handler(ncstat)
316
317!-----------------------------------------------------------------------
318!
319!     write grid data
320!
321!-----------------------------------------------------------------------
322
323      ncstat = nf_put_var_int(nc_grid_id, nc_griddims_id, grid_dims)
324      call netcdf_error_handler(ncstat)
325
326      ncstat = nf_put_var_int(nc_grid_id, nc_grdimask_id, grid_imask)
327      call netcdf_error_handler(ncstat)
328
329      ncstat = nf_put_var_double(nc_grid_id, nc_grdcntrlat_id, 
330     &                           grid_center_lat)
331      call netcdf_error_handler(ncstat)
332
333      ncstat = nf_put_var_double(nc_grid_id, nc_grdcntrlon_id, 
334     &                           grid_center_lon)
335      call netcdf_error_handler(ncstat)
336
337      ncstat = nf_put_var_double(nc_grid_id, nc_grdcrnrlat_id, 
338     &                           grid_corner_lat)
339      call netcdf_error_handler(ncstat)
340
341      ncstat = nf_put_var_double(nc_grid_id, nc_grdcrnrlon_id, 
342     &                           grid_corner_lon)
343      call netcdf_error_handler(ncstat)
344
345      ncstat = nf_close(nc_grid_id)
346      call netcdf_error_handler(ncstat)
347
348!-----------------------------------------------------------------------
349
350      end program convert_gauss
351
352!***********************************************************************
353
354      subroutine gquad(l,root,w)
355
356!-----------------------------------------------------------------------
357!
358!     This subroutine finds the l roots (in theta) and gaussian weights
359!     associated with the legendre polynomial of degree l > 1.
360!
361!-----------------------------------------------------------------------
362
363      use kinds_mod
364      use constants
365
366      implicit none
367
368!-----------------------------------------------------------------------
369!
370!     intent(in)
371!
372!-----------------------------------------------------------------------
373
374      integer (kind=int_kind), intent(in) :: l
375
376!-----------------------------------------------------------------------
377!
378!     intent(out)
379!
380!-----------------------------------------------------------------------
381
382      real (kind=dbl_kind), dimension(l), intent(out) ::
383     &      root, w
384
385!-----------------------------------------------------------------------
386!
387!     local
388!
389!-----------------------------------------------------------------------
390
391      integer (kind=int_kind) :: l1, l2, l22, l3, k, i, j
392
393      real (kind=dbl_kind) :: 
394     &     del,co,p1,p2,p3,t1,t2,slope,s,c,pp1,pp2,p00
395
396!-----------------------------------------------------------------------
397!
398!     Define useful constants.
399!
400!-----------------------------------------------------------------------
401
402      del= pi/float(4*l)
403      l1 = l+1
404      co = float(2*l+3)/float(l1**2)
405      p2 = 1.0
406      t2 = -del
407      l2 = l/2
408      k = 1
409      p00 = one/sqrt(two)
410
411!-----------------------------------------------------------------------
412!
413!     Start search for each root by looking for crossing point.
414!
415!-----------------------------------------------------------------------
416
417      do i=1,l2
418   10    t1 = t2
419         t2 = t1+del
420         p1 = p2
421         s = sin(t2)
422         c = cos(t2)
423         pp1 = 1.0
424         p3 = p00
425         do j=1,l1
426            pp2 = pp1
427            pp1 = p3
428            p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-
429     &           sqrt(float((2*j+1)*(j-1)*(j-1))/
430     &           float((2*j-3)*j*j))*pp2
431         end do
432         p2 = pp1
433         if ((k*p2).gt.0) goto 10
434
435!-----------------------------------------------------------------------
436!
437!        Now converge using Newton-Raphson.
438!
439!-----------------------------------------------------------------------
440
441         k = -k
442   20    continue
443            slope = (t2-t1)/(p2-p1)
444            t1 = t2
445            t2 = t2-slope*p2
446            p1 = p2
447            s = sin(t2)
448            c = cos(t2)
449            pp1 = 1.0
450            p3 = p00
451            do j=1,l1
452               pp2 = pp1
453               pp1 = p3
454               p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-
455     &              sqrt(float((2*j+1)*(j-1)*(j-1))/
456     &              float((2*j-3)*j*j))*pp2
457            end do
458            p2 = pp1
459         if (abs(p2).gt.1.e-10) goto 20
460         root(i) = t2
461         w(i) = co*(sin(t2)/p3)**2
462      end do
463
464!-----------------------------------------------------------------------
465!
466!     If l is odd, take care of odd point.
467!
468!-----------------------------------------------------------------------
469
470      l22 = 2*l2
471      if (l22 .ne. l) then
472         l2 = l2+1
473         t2 = pi/2.0
474         root(l2) = t2
475         s = sin(t2)
476         c = cos(t2)
477         pp1 = 1.0
478         p3 = p00
479         do j=1,l1
480            pp2 = pp1
481            pp1 = p3
482            p3 = 2.0*sqrt((float(j**2)-0.250)/float(j**2))*c*pp1-
483     &           sqrt(float((2*j+1)*(j-1)*(j-1))/
484     &           float((2*j-3)*j*j))*pp2
485         end do
486         p2 = pp1
487         w(l2) = co/p3**2
488      endif
489
490!-----------------------------------------------------------------------
491!
492!     Use symmetry to compute remaining roots and weights.
493!
494!-----------------------------------------------------------------------
495
496      l3 = l2+1
497      do i=l3,l
498         root(i) = pi-root(l-i+1)
499         w(i) = w(l-i+1)
500      end do
501
502!-----------------------------------------------------------------------
503
504      end subroutine gquad
505
506!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.