source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/lib/mct/testsystem/testall/convertPOPT.F90 @ 6331

Last change on this file since 6331 was 6331, checked in by aclsce, 17 months ago

Moved oasis-mct_5.0 in oasis3-mct/branches directory.

File size: 14.6 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.F90,v 1.9 2004-06-02 23:25:50 eong Exp $
9!     CVS $Name:  $
10!
11!     Copyright (c) 1997, 1998 the Regents of the University of
12!       California.
13!
14!     Unless otherwise indicated, this software has been authored
15!     by an employee or employees of the University of California,
16!     operator of the Los Alamos National Laboratory under Contract
17!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
18!     Government has rights to use, reproduce, and distribute this
19!     software.  The public may copy and use this software without
20!     charge, provided that this Notice and any statement of authorship
21!     are reproduced on all copies.  Neither the Government nor the
22!     University makes any warranty, express or implied, or assumes
23!     any liability or responsibility for the use of this software.
24!
25!***********************************************************************
26
27      subroutine convertPOPT(GGrid, grid_file_in, grid_topo_in, nx, ny)
28
29!-----------------------------------------------------------------------
30!
31!     This file converts a POP grid.dat file to a remapping grid file.
32!
33!-----------------------------------------------------------------------
34
35      use m_AttrVect,only    : AttrVect
36      use m_GeneralGrid,only : MCT_GGrid_init => init
37      use m_GeneralGrid,only : MCT_GGrid_indexIA => indexIA
38      use m_GeneralGrid,only : MCT_GGrid_indexRA => indexRA
39      use m_GeneralGrid,only : GeneralGrid
40      use m_stdio
41      use m_ioutil
42      use m_die
43
44
45      implicit none
46
47!-----------------------------------------------------------------------
48!
49!     variables that describe the grid
50!       4/3       nx = 192, ny = 128
51!       2/3 (mod) nx = 384, ny = 288
52!       x3p Greenland DP nx = 100, ny = 116
53!       x2p Greenland DP nx = 160, ny = 192
54!       x1p Greenland DP nx = 320, ny = 384
55!
56!-----------------------------------------------------------------------
57
58      type(GeneralGrid),       intent(out) :: GGrid
59      character (len=*),       intent(in)  :: grid_file_in
60      character (len=*),       intent(in)  :: grid_topo_in
61      integer,                 intent(in)  :: nx
62      integer,                 intent(in)  :: ny
63
64      integer :: grid_size
65
66      integer, parameter :: &
67                   grid_rank = 2, &
68                   grid_corners = 4
69
70      integer, dimension(2) :: &
71                   grid_dims   ! size of each dimension
72
73!-----------------------------------------------------------------------
74!
75!     grid coordinates and masks
76!
77!-----------------------------------------------------------------------
78
79!:: NOTE: The following kind specifiers are needed to read the proper
80!:: values for the POP grid files. The subsequent type conversions
81!:: on these variables may pose a risk.
82
83      integer(kind(1)), dimension(:), allocatable :: &
84                   grid_imask
85
86      real, dimension(:), allocatable :: &
87                   grid_area      ,  &! area as computed in POP
88                   grid_center_lat,  &! lat/lon coordinates for
89                   grid_center_lon   ! each grid center in radians
90
91      real(selected_real_kind(13)), dimension(:,:), allocatable :: &
92                   grid_corner_lat,  &! lat/lon coordinates for
93                   grid_corner_lon   ! each grid corner in radians
94
95      real(selected_real_kind(13)), dimension(:,:), allocatable :: &
96                   HTN, HTE          ! T-cell grid lengths
97
98!-----------------------------------------------------------------------
99!
100!     defined constants
101!
102!-----------------------------------------------------------------------
103
104      real(selected_real_kind(13)), parameter ::  &
105            zero   = 0.0,           &
106            one    = 1.0,           &
107            two    = 2.0,           &
108            three  = 3.0,           &
109            four   = 4.0,           &
110            five   = 5.0,           &
111            half   = 0.5,           &
112            quart  = 0.25,          &
113            bignum = 1.e+20,        &
114            tiny   = 1.e-14,        &
115            pi     = 3.14159265359, &
116            pi2    = two*pi,        &
117            pih    = half*pi
118
119      real(selected_real_kind(13)), parameter ::  &
120                   radius    = 6.37122e8 ,        &  ! radius of Earth (cm)
121                   area_norm = one/(radius*radius)
122
123!-----------------------------------------------------------------------
124!
125!     other local variables
126!
127!-----------------------------------------------------------------------
128
129      character(len=*),parameter :: myname_= 'convertPOPT'
130
131      integer :: i, j, k, n, p, q, r, ier
132
133      integer :: iunit, ocn_add, im1, jm1, np1, np2
134
135      integer :: center_lat, center_lon, &
136                 corner_lat, corner_lon, &
137                 imask, area
138
139      real :: tmplon, dlat, dxt, dyt
140
141      real :: x1, x2, x3, x4, &
142              y1, y2, y3, y4, &
143              z1, z2, z3, z4, &
144              tx, ty, tz, da
145
146      grid_size = nx*ny
147
148      allocate(grid_imask(grid_size), &
149               grid_area(grid_size), &
150               grid_center_lat(grid_size), &
151               grid_center_lon(grid_size), &
152               grid_corner_lat(grid_corners,grid_size), &
153               grid_corner_lon(grid_corners,grid_size), &
154               HTN(nx,ny), &
155               HTE(nx,ny), &
156               stat=ier)
157
158      if(ier/=0) call die(myname_,"allocate(grid_imask... ", ier)
159
160!-----------------------------------------------------------------------
161!
162!     read in grid info
163!     lat/lon info is on velocity points which correspond
164!     to the NE corner (in logical space) of the grid cell.
165!
166!-----------------------------------------------------------------------
167
168      iunit = luavail()
169
170      open(unit=iunit, file=trim(grid_topo_in), status='old', &
171           form='unformatted', access='direct', recl=grid_size*4)
172
173      read (unit=iunit,rec=1) grid_imask
174
175      call luflush(iunit)
176
177      iunit = luavail()
178#if SYSSUPERUX || SYSOSF1
179      open(unit=iunit, file=trim(grid_file_in), status='old', &
180           form='unformatted', access='direct', recl=grid_size*2)
181#else
182      open(unit=iunit, file=trim(grid_file_in), status='old', &
183           form='unformatted', access='direct', recl=grid_size*8)
184#endif
185
186      read (unit=iunit, rec=1) grid_corner_lat(3,:)
187      read (unit=iunit, rec=2) grid_corner_lon(3,:)
188      read (unit=iunit, rec=3) HTN
189      read (unit=iunit, rec=4) HTE
190      call luflush(iunit)
191
192!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
193!::::::::::::TEST DIAGNOSTICS::::::::::::::::::::::::::::::::::
194!::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
195
196      k=0
197      do j=1,grid_size
198         if(grid_imask(j)==0) k=k+1
199      enddo
200
201      write(stdout,*) "CONVERTPOPT: NUM_ZEROES(GRID_IMASK), SUM(GRID_IMASK)",&
202           k, sum(grid_imask)
203
204     write(stdout,*) "CONVERTPOPT: GRID_CORNER_LAT VALUES = ", &
205          grid_corner_lat(3,1:10)
206
207     write(stdout,*) "CONVERTPOPT: GRID_CORNER_LON VALUES = ", &
208          grid_corner_lon(3,1:10)
209
210     write(stdout,*) "CONVERTPOPT: HTN VALUES = ", &
211          HTN(1,1:10)
212
213     write(stdout,*) "CONVERTPOPT: HTE VALUES = ", &
214          HTE(1,1:10)
215
216!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
217
218      grid_dims(1) = nx
219      grid_dims(2) = ny
220
221!-----------------------------------------------------------------------
222!
223!     convert KMT field to integer grid mask
224!
225!-----------------------------------------------------------------------
226
227      grid_imask = min(grid_imask, 1)
228
229!-----------------------------------------------------------------------
230!
231!     compute remaining corners
232!
233!-----------------------------------------------------------------------
234
235      do j=1,ny
236        do i=1,nx
237          ocn_add = (j-1)*nx + i
238          if (i .ne. 1) then
239            im1 = ocn_add - 1
240          else
241            im1 = ocn_add + nx - 1
242          endif
243
244          grid_corner_lat(4,ocn_add) = grid_corner_lat(3,im1)
245          grid_corner_lon(4,ocn_add) = grid_corner_lon(3,im1)
246        end do
247      end do
248
249      do j=2,ny
250        do i=1,nx
251          ocn_add = (j-1)*nx + i
252          jm1 = (j-2)*nx + i
253
254          grid_corner_lat(2,ocn_add) = grid_corner_lat(3,jm1)
255          grid_corner_lat(1,ocn_add) = grid_corner_lat(4,jm1)
256
257          grid_corner_lon(2,ocn_add) = grid_corner_lon(3,jm1)
258          grid_corner_lon(1,ocn_add) = grid_corner_lon(4,jm1)
259        end do
260      end do
261
262!-----------------------------------------------------------------------
263!
264!     mock up the lower row boundaries
265!
266!-----------------------------------------------------------------------
267
268      do i=1,nx
269        dlat = grid_corner_lat(1,i+2*nx) - grid_corner_lat(1,i+nx)
270        grid_corner_lat(1,i) = grid_corner_lat(1,i+nx) - dlat
271        grid_corner_lat(1,i) = max(grid_corner_lat(1,i), -pih + tiny)
272
273        dlat = grid_corner_lat(2,i+2*nx) - grid_corner_lat(2,i+nx)
274        grid_corner_lat(2,i) = grid_corner_lat(2,i+nx) - dlat
275        grid_corner_lat(2,i) = max(grid_corner_lat(2,i), -pih + tiny)
276
277        grid_corner_lon(1,i) = grid_corner_lon(4,i)
278        grid_corner_lon(2,i) = grid_corner_lon(3,i)
279      end do
280
281!-----------------------------------------------------------------------
282!
283!     correct for 0,2pi longitude crossings
284!
285!-----------------------------------------------------------------------
286
287      do ocn_add=1,grid_size
288        if (grid_corner_lon(1,ocn_add) > pi2) &
289            grid_corner_lon(1,ocn_add) = &
290            grid_corner_lon(1,ocn_add) - pi2
291        if (grid_corner_lon(1,ocn_add) < 0.0) &
292            grid_corner_lon(1,ocn_add) = &
293            grid_corner_lon(1,ocn_add) + pi2
294        do n=2,grid_corners
295          tmplon = grid_corner_lon(n  ,ocn_add) -  &
296                   grid_corner_lon(n-1,ocn_add)
297          if (tmplon < -three*pih) grid_corner_lon(n,ocn_add) = &
298                                   grid_corner_lon(n,ocn_add) + pi2
299          if (tmplon >  three*pih) grid_corner_lon(n,ocn_add) = &
300                                   grid_corner_lon(n,ocn_add) - pi2
301        end do
302      end do
303
304!-----------------------------------------------------------------------
305!
306!     compute ocean cell centers by averaging corner values
307!
308!-----------------------------------------------------------------------
309
310      do ocn_add=1,grid_size
311         z1 = cos(grid_corner_lat(1,ocn_add))
312         x1 = cos(grid_corner_lon(1,ocn_add))*z1
313         y1 = sin(grid_corner_lon(1,ocn_add))*z1
314         z1 = sin(grid_corner_lat(1,ocn_add))
315
316         z2 = cos(grid_corner_lat(2,ocn_add))
317         x2 = cos(grid_corner_lon(2,ocn_add))*z2
318         y2 = sin(grid_corner_lon(2,ocn_add))*z2
319         z2 = sin(grid_corner_lat(2,ocn_add))
320
321         z3 = cos(grid_corner_lat(3,ocn_add))
322         x3 = cos(grid_corner_lon(3,ocn_add))*z3
323         y3 = sin(grid_corner_lon(3,ocn_add))*z3
324         z3 = sin(grid_corner_lat(3,ocn_add))
325
326         z4 = cos(grid_corner_lat(4,ocn_add))
327         x4 = cos(grid_corner_lon(4,ocn_add))*z4
328         y4 = sin(grid_corner_lon(4,ocn_add))*z4
329         z4 = sin(grid_corner_lat(4,ocn_add))
330
331         tx = (x1+x2+x3+x4)/4.0
332         ty = (y1+y2+y3+y4)/4.0
333         tz = (z1+z2+z3+z4)/4.0
334         da = sqrt(tx**2+ty**2+tz**2)
335
336         tz = tz/da
337                                ! grid_center_lon in radians
338         grid_center_lon(ocn_add) = 0.0
339         if (tx .ne. 0.0 .or. ty .ne. 0.0) &
340              grid_center_lon(ocn_add) = atan2(ty,tx)
341                                ! grid_center_lat in radians
342         grid_center_lat(ocn_add) = asin(tz)
343
344      end do
345
346      ! j=1: linear approximation
347      n = 0
348      do i=1,nx
349         n   = n + 1
350         np1 = n + nx
351         np2 = n + 2*nx
352         grid_center_lon(n) = grid_center_lon(np1)
353         grid_center_lat(n) = 2.0*grid_center_lat(np1) - &
354              grid_center_lat(np2)
355      end do
356
357      do ocn_add=1,grid_size
358        if (grid_center_lon(ocn_add) > pi2) &
359            grid_center_lon(ocn_add) = grid_center_lon(ocn_add) - pi2
360        if (grid_center_lon(ocn_add) < 0.0) &
361            grid_center_lon(ocn_add) = grid_center_lon(ocn_add) + pi2
362      enddo
363
364!-----------------------------------------------------------------------
365!
366!     compute cell areas in same way as POP
367!
368!-----------------------------------------------------------------------
369
370      n = 0
371      do j=1,ny
372        if (j > 1) then
373          jm1 = j-1
374        else
375          jm1 = 1
376        endif
377        do i=1,nx
378          if (i > 1) then
379            im1 = i-1
380          else
381            im1 = nx
382          endif
383
384          n = n+1
385
386          dxt = half*(HTN(i,j) + HTN(i,jm1))
387          dyt = half*(HTE(i,j) + HTE(im1,j))
388          if (dxt == zero) dxt=one
389          if (dyt == zero) dyt=one
390
391          grid_area(n) = dxt*dyt*area_norm
392        end do
393      end do
394
395!-----------------------------------------------------------------------
396!
397!     intialize GeneralGrid
398!
399!-----------------------------------------------------------------------
400
401      call MCT_GGrid_init(GGrid=GGrid, &
402                          CoordChars="grid_center_lat:&
403                                     &grid_center_lon", &
404                          WeightChars="grid_area", &
405                          OtherChars="grid_corner_lat_1:&
406                                     &grid_corner_lat_2:&
407                                     &grid_corner_lat_3:&
408                                     &grid_corner_lat_4:&
409                                     &grid_corner_lon_1:&
410                                     &grid_corner_lon_2:&
411                                     &grid_corner_lon_3:&
412                                     &grid_corner_lon_4", &
413                          IndexChars="grid_imask", &
414                          lsize=grid_size)
415
416      center_lat = MCT_GGrid_indexRA(GGrid,'grid_center_lat')
417      center_lon = MCT_GGrid_indexRA(GGrid,'grid_center_lon')
418      corner_lat = MCT_GGrid_indexRA(GGrid,'grid_corner_lat_1')
419      corner_lon = MCT_GGrid_indexRA(GGrid,'grid_corner_lon_1')
420      area       = MCT_GGrid_indexRA(GGrid,'grid_area')
421      imask      = MCT_GGrid_indexIA(GGrid,'grid_imask')
422
423      GGrid%data%rattr(center_lat,1:grid_size) = &
424           grid_center_lat(1:grid_size)
425      GGrid%data%rattr(center_lon,1:grid_size) = &
426           grid_center_lon(1:grid_size)
427      GGrid%data%rattr(area,1:grid_size) = &
428           grid_area(1:grid_size)
429      GGrid%data%iattr(imask,1:grid_size) = &
430           grid_imask(1:grid_size)
431
432      do p = 1,grid_corners
433         GGrid%data%rattr(corner_lat+p-1,1:grid_size) = &
434              grid_corner_lat(p,1:grid_size)
435         GGrid%data%rattr(corner_lon+p-1,1:grid_size) = &
436           grid_corner_lon(p,1:grid_size)
437      enddo
438
439      deallocate(grid_imask, grid_area,            &
440                 grid_center_lat, grid_center_lon, &
441                 grid_corner_lat, grid_corner_lon, &
442                 HTN, HTE, stat=ier)
443
444      if(ier/=0) call die(myname_,"deallocate(grid_imask... ", ier)
445
446
447!***********************************************************************
448
449    end subroutine convertPOPT
450
451!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
452
453
454
Note: See TracBrowser for help on using the repository browser.