source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/lib/psmile/src/mod_oasis_grid.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: 57.8 KB
Line 
1
2!> OASIS grid data and methods
3
4!> These interfaces support both grid data specified globally on the root task
5!> as required in Oasis3 and grid data decomposed across tasks.  If grid data
6!> is decomposed across tasks, the optional partid argument must be specified
7!> when it exists in the interface.
8
9MODULE mod_oasis_grid
10!-----------------------------------------------------------------------
11! BOP
12!
13! !MODULE:  mod_prism_grid
14! !REMARKS:
15!
16! **********************
17! THIS SHOULD BE CALLED BY A SINGLE PE ACCORDING TO THE OASIS3
18! STANDARD.  THE DATA IS GLOBAL.
19! **********************
20!
21! !REVISION HISTORY:
22!
23!
24! !PUBLIC MEMBER FUNCTIONS:
25!
26!      subroutine oasis_start_grids_writing(iwrite)
27!             This subroutine initializes grid writing by receiving a
28!             starting command from OASIS.
29!
30!      subroutine oasis_write_grid(cgrid, nx, ny, lon, lat, part_id)
31!             This subroutine writes longitudes and latitudes for a model
32!             grid.  part_id is optional and indicates decomposed data is provided
33!
34!      subroutine oasis_write_corner(cgrid, nx, ny, nc, clon, clat, part_id)
35!             This subroutine writes the longitudes and latitudes of the
36!             grid cell corners.  part_id is optional and indicates decomposed data
37!             is provided
38!
39!      subroutine oasis_write_mask(cgrid, nx, ny, mask, part_id, companion)
40!             This subroutine writes the mask for a model grid.  part_id is optional
41!             and indicates decomposed data is provided
42!
43!      subroutine oasis_write_area(cgrid, nx, ny, area, part_id)
44!             This subroutine writes the grid cell areas for a model grid.  part_id
45!             is optional and indicates decomposed data is provided
46!
47!      subroutine oasis_write_frac(cgrid, nx, ny, frac, part_id, companion)
48!             This subroutine writes the grid cell fracs for a model grid.  part_id
49!             is optional and indicates decomposed data is provided
50!
51!      subroutine oasis_terminate_grids_writing()
52!             This subroutine terminates grid writing by sending a flag
53!             to OASIS, stating the all needed grid information was written.
54!       
55
56! !USES:
57  USE mod_oasis_data
58  USE mod_oasis_io
59  USE mod_oasis_sys
60  USE mod_oasis_part
61  USE mod_oasis_mpi, only: oasis_mpi_min, oasis_mpi_max, oasis_mpi_bcast, oasis_mpi_barrier, &
62                           oasis_mpi_chkerr, oasis_mpi_reducelists
63  USE mod_oasis_timer
64  USE mct_mod
65 
66  implicit none
67
68  private
69
70  public oasis_start_grids_writing
71  public oasis_write_grid
72  public oasis_write_angle
73  public oasis_write_corner
74  public oasis_write_mask
75  public oasis_write_area
76  public oasis_write_frac
77  public oasis_terminate_grids_writing 
78  public oasis_write2files
79  public oasis_print_grid_data
80
81#include "oasis_os.h"
82
83  !> Generic interface to support writing 4 or 8 byte reals
84  interface oasis_write_grid
85#ifndef __NO_4BYTE_REALS
86     module procedure oasis_write_grid_r4
87#endif
88     module procedure oasis_write_grid_r8
89  end interface
90
91  !> Generic interface to support writing 4 or 8 byte reals
92  interface oasis_write_angle
93#ifndef __NO_4BYTE_REALS
94     module procedure oasis_write_angle_r4
95#endif
96     module procedure oasis_write_angle_r8
97  end interface
98
99  !> Generic interface to support writing 4 or 8 byte reals
100  interface oasis_write_corner
101#ifndef __NO_4BYTE_REALS
102     module procedure oasis_write_corner_r4
103#endif
104     module procedure oasis_write_corner_r8
105  end interface
106
107  !> Generic interface to support writing 4 or 8 byte reals
108  interface oasis_write_area
109#ifndef __NO_4BYTE_REALS
110     module procedure oasis_write_area_r4
111#endif
112     module procedure oasis_write_area_r8
113  end interface
114
115  !> Generic interface to support writing 4 or 8 byte reals
116  interface oasis_write_frac
117#ifndef __NO_4BYTE_REALS
118     module procedure oasis_write_frac_r4
119#endif
120     module procedure oasis_write_frac_r8
121  end interface
122
123  !--- datatypes ---
124  public :: prism_grid_type    !< Grid datatype
125
126  integer(kind=ip_intwp_p),parameter :: mgrid = 100  !< maximum number of grids allowed
127  integer(kind=ip_intwp_p),save :: writing_grids_call=0
128
129  !> Model grid data for creating mapping data and conserving fields
130  type prism_grid_type
131     character(len=ic_med)  :: gridname   !< grid name
132     integer(kind=ip_i4_p)  :: partid     !< partition ID
133     integer(kind=ip_i4_p)  :: nx         !< global nx size
134     integer(kind=ip_i4_p)  :: ny         !< global ny size
135     integer(kind=ip_i4_p)  :: nc         !< number of corners per gridcell
136     logical                :: grid_set   !< flag to track user calls for grid
137     logical                :: corner_set !< flag to track user calls for corner
138     logical                :: angle_set  !< flag to track user calls for angle
139     logical                :: area_set   !< flag to track user calls for area
140     logical                :: frac_set   !< flag to track user calls for frac
141     logical                :: mask_set   !< flag to track user calls for mask
142     logical                :: written    !< flag to indicate grid has been written
143     logical                :: terminated !< flag to indicate user grid calls complete
144     character(len=ic_med)  :: mask_companion   !< mask companion grid name
145     character(len=ic_med)  :: frac_companion   !< frac companion grid name
146     real(kind=ip_realwp_p),allocatable :: lon(:,:)     !< user specified longitudes
147     real(kind=ip_realwp_p),allocatable :: lat(:,:)     !< user specified latitudes
148     real(kind=ip_realwp_p),allocatable :: clon(:,:,:)  !< user specified corner longitudes
149     real(kind=ip_realwp_p),allocatable :: clat(:,:,:)  !< user specified corner latitudes
150     real(kind=ip_realwp_p),allocatable :: angle(:,:)   !< user specified angle
151     real(kind=ip_realwp_p),allocatable :: area(:,:)    !< user specified area
152     real(kind=ip_realwp_p),allocatable :: frac(:,:)    !< user specified frac
153     integer(kind=ip_i4_p) ,allocatable :: mask(:,:)    !< user specified mask
154  end type prism_grid_type
155
156  integer(kind=ip_intwp_p),public,save :: prism_ngrid = 0  !< counter for grids
157  type(prism_grid_type),public,save :: prism_grid(mgrid)   !< array of grid datatypes
158  logical, parameter :: local_timers_on = .false.
159
160#ifdef use_netCDF
161#include <netcdf.inc>
162#endif
163
164!---------------------------------------------------------------------------
165
166CONTAINS
167
168!--------------------------------------------------------------------------
169
170!> Print grid information to log file.
171
172    SUBROUTINE oasis_print_grid_data()
173
174    implicit none
175 
176    !-------------------------------------------------
177    integer(kind=ip_intwp_p) :: n
178    character(len=*),parameter :: subname = '(oasis_print_grid_data)'
179    !-------------------------------------------------
180
181    call oasis_debug_enter(subname)
182    if (local_timers_on) call oasis_timer_start(trim(subname))
183
184    if (OASIS_debug >= 15) then
185    do n = 1,prism_ngrid
186       write(nulprt,*) ' '
187       write(nulprt,*) subname,trim(prism_grid(n)%gridname),' size',prism_grid(n)%nx,prism_grid(n)%ny
188       write(nulprt,*) subname,trim(prism_grid(n)%gridname),' set ',prism_grid(n)%grid_set, &
189                       prism_grid(n)%corner_set, prism_grid(n)%angle_set, prism_grid(n)%area_set, &
190                       prism_grid(n)%frac_set, prism_grid(n)%mask_set, &
191                       trim(prism_grid(n)%mask_companion), ' ', trim(prism_grid(n)%frac_companion)
192       if (prism_grid(n)%partid > 0 .and. prism_grid(n)%partid < prism_npart) then
193          write(nulprt,*) subname,'partid ',trim(prism_grid(n)%gridname),prism_grid(n)%partid, &
194                       trim(prism_part(prism_grid(n)%partid)%partname)
195       else
196          write(nulprt,*) subname,'partid ',trim(prism_grid(n)%gridname),prism_grid(n)%partid
197       endif
198       if (prism_grid(n)%grid_set) write(nulprt,*) subname,trim(prism_grid(n)%gridname),' lon  ', &
199                                   size(prism_grid(n)%lon,dim=1),size(prism_grid(n)%lon,dim=2), &
200                                   minval(prism_grid(n)%lon),maxval(prism_grid(n)%lon)
201       if (prism_grid(n)%grid_set) write(nulprt,*) subname,trim(prism_grid(n)%gridname),' lat  ', &
202                                   size(prism_grid(n)%lat,dim=1),size(prism_grid(n)%lat,dim=2), &
203                                   minval(prism_grid(n)%lat),maxval(prism_grid(n)%lat)
204       if (prism_grid(n)%corner_set) write(nulprt,*) subname,trim(prism_grid(n)%gridname),' clon ', &
205                                   size(prism_grid(n)%clon,dim=1),size(prism_grid(n)%clon,dim=2), &
206                                   minval(prism_grid(n)%clon),maxval(prism_grid(n)%clon)
207       if (prism_grid(n)%corner_set) write(nulprt,*) subname,trim(prism_grid(n)%gridname),' clat ', &
208                                   size(prism_grid(n)%clat,dim=1),size(prism_grid(n)%clat,dim=2), &
209                                   minval(prism_grid(n)%clat),maxval(prism_grid(n)%clat)
210       if (prism_grid(n)%angle_set) write(nulprt,*) subname,trim(prism_grid(n)%gridname),' angl ', &
211                                   size(prism_grid(n)%angle,dim=1),size(prism_grid(n)%angle,dim=2), &
212                                   minval(prism_grid(n)%angle),maxval(prism_grid(n)%angle)
213       if (prism_grid(n)%mask_set) write(nulprt,*) subname,trim(prism_grid(n)%gridname),' mask ', &
214                                   size(prism_grid(n)%mask,dim=1),size(prism_grid(n)%mask,dim=2), &
215                                   minval(prism_grid(n)%mask),maxval(prism_grid(n)%mask)
216       if (prism_grid(n)%area_set) write(nulprt,*) subname,trim(prism_grid(n)%gridname),' area ', &
217                                   size(prism_grid(n)%area,dim=1),size(prism_grid(n)%area,dim=2), &
218                                   minval(prism_grid(n)%area),maxval(prism_grid(n)%area)
219       if (prism_grid(n)%frac_set) write(nulprt,*) subname,trim(prism_grid(n)%gridname),' frac ', &
220                                   size(prism_grid(n)%frac,dim=1),size(prism_grid(n)%frac,dim=2), &
221                                   minval(prism_grid(n)%frac),maxval(prism_grid(n)%frac)
222    enddo
223    endif
224
225    if (local_timers_on) call oasis_timer_stop(trim(subname))
226    call oasis_debug_exit(subname)
227
228  END SUBROUTINE oasis_print_grid_data
229
230!--------------------------------------------------------------------------
231
232!> User interface to initialize grid writing
233
234    SUBROUTINE oasis_start_grids_writing(iwrite)
235
236    implicit none
237 
238    integer(kind=ip_intwp_p), intent (OUT) :: iwrite !< flag, obsolete
239
240    !-------------------------------------------------
241    character(len=*),parameter :: subname = '(oasis_start_grids_writing)'
242    !-------------------------------------------------
243
244    call oasis_debug_enter(subname)
245    if (local_timers_on) call oasis_timer_start(trim(subname))
246
247    if (OASIS_debug >= 15) then
248       write(nulprt,*) subname,' prism_ngrid = ',prism_ngrid
249    endif
250
251    if (prism_ngrid == 0) then  ! first call
252       prism_grid(:)%gridname   = 'unSet'
253       prism_grid(:)%nx         = -1
254       prism_grid(:)%ny         = -1
255       prism_grid(:)%grid_set   = .false.
256       prism_grid(:)%corner_set = .false.
257       prism_grid(:)%angle_set  = .false.
258       prism_grid(:)%area_set   = .false.
259       prism_grid(:)%frac_set   = .false.
260       prism_grid(:)%mask_set   = .false.
261       prism_grid(:)%written    = .false.
262       prism_grid(:)%terminated = .false.
263       prism_grid(:)%partid     = -1
264       prism_grid(:)%mask_companion = 'undefined'
265       prism_grid(:)%frac_companion = 'undefined'
266    endif
267    iwrite = 1   ! just set grids are needed always
268    writing_grids_call=1
269
270    if (local_timers_on) call oasis_timer_stop(trim(subname))
271    call oasis_debug_exit(subname)
272
273  END SUBROUTINE oasis_start_grids_writing
274
275!--------------------------------------------------------------------------
276
277!> User interface to set latitudes and longitudes for 8 byte reals
278
279    SUBROUTINE oasis_write_grid_r8(cgrid, nx, ny, lon, lat, partid)
280
281    !-------------------------------------------------
282    ! Routine to create a new grids file or to add a grid description to an
283    ! existing grids file.
284    !-------------------------------------------------
285
286    implicit none
287
288    character(len=*),         intent (in) :: cgrid      !< grid name
289    integer(kind=ip_intwp_p), intent (in) :: nx         !< global nx size
290    integer(kind=ip_intwp_p), intent (in) :: ny         !< global ny size
291    real(kind=ip_double_p),   intent (in) :: lon(:,:)   !< longitudes
292    real(kind=ip_double_p),   intent (in) :: lat(:,:)   !< latitudes
293    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
294    !-------------------------------------------------
295    integer(kind=ip_intwp_p) :: GRIDID
296    integer(kind=ip_intwp_p) :: ierror
297    integer(kind=ip_intwp_p) :: lnx,lny
298    character(len=*),parameter :: subname = '(oasis_write_grid_r8)'
299    !-------------------------------------------------
300
301    call oasis_debug_enter(subname)
302    if (local_timers_on) call oasis_timer_start(trim(subname))
303
304    if (OASIS_debug >= 15) then
305       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
306    endif
307
308    call oasis_findgrid(cgrid,nx,ny,gridID)
309
310    lnx = size(lon,dim=1)
311    lny = size(lon,dim=2)
312
313    allocate(prism_grid(gridID)%lon(lnx,lny),stat=ierror)
314    IF (ierror /= 0) WRITE(nulprt,*) subname,' model :',compid,' proc :',&
315                                     mpi_rank_local,' WARNING lon alloc'
316
317    lnx = size(lat,dim=1)
318    lny = size(lat,dim=2)
319
320    allocate(prism_grid(gridID)%lat(lnx,lny),stat=ierror)
321    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
322                                     mpi_rank_local,' WARNING lat alloc'
323
324    prism_grid(gridID)%lon = lon
325    prism_grid(gridID)%lat = lat
326    prism_grid(gridID)%grid_set = .true.
327
328    if (present(partid)) then
329       if (prism_grid(gridID)%partid > 0 .and. prism_grid(gridID)%partid /= partid) then
330          write(nulprt,*) subname,estr,'partid inconsistency',gridID,prism_grid(gridID)%partid,partid
331          call oasis_abort(file=__FILE__,line=__LINE__)
332       endif
333       prism_grid(gridID)%partid = partid
334       if (OASIS_debug >= 15) then
335          write(nulprt,*) subname,' partid = ',trim(cgrid),partid
336       endif
337    endif
338
339    if (local_timers_on) call oasis_timer_stop(trim(subname))
340    call oasis_debug_exit(subname)
341
342  END SUBROUTINE oasis_write_grid_r8
343
344!--------------------------------------------------------------------------
345
346!> User interface to set latitudes and longitudes for 4 byte reals
347
348    SUBROUTINE oasis_write_grid_r4(cgrid, nx, ny, lon, lat, partid)
349
350    !-------------------------------------------------
351    ! Routine to create a new grids file or to add a grid description to an
352    ! existing grids file.
353    !-------------------------------------------------
354
355    implicit none
356
357    character(len=*),         intent (in) :: cgrid      !< grid name
358    integer(kind=ip_intwp_p), intent (in) :: nx         !< global nx size
359    integer(kind=ip_intwp_p), intent (in) :: ny         !< global ny size
360    real(kind=ip_single_p),   intent (in) :: lon(:,:)   !< longitudes
361    real(kind=ip_single_p),   intent (in) :: lat(:,:)   !< latitudes
362    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
363    !-------------------------------------------------
364    real(kind=ip_double_p), allocatable :: lon8(:,:)
365    real(kind=ip_double_p), allocatable :: lat8(:,:)
366    integer(kind=ip_intwp_p) :: ierror
367    integer(kind=ip_intwp_p) :: lpartid
368    integer(kind=ip_intwp_p) :: lnx,lny
369    character(len=*),parameter :: subname = '(oasis_write_grid_r4)'
370    !-------------------------------------------------
371
372    call oasis_debug_enter(subname)
373    if (local_timers_on) call oasis_timer_start(trim(subname))
374
375    if (OASIS_debug >= 15) then
376       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
377    endif
378
379    lpartid = -1
380    if (present(partid)) then
381       lpartid = partid
382    endif
383    if (OASIS_debug >= 15) then
384       write(nulprt,*) subname,' partid = ',trim(cgrid),lpartid
385    endif
386
387    lnx = size(lon,dim=1)
388    lny = size(lon,dim=2)
389
390    allocate(lon8(lnx,lny),stat=ierror)
391    IF (ierror /= 0) WRITE(nulprt,*) subname,' model :',compid,' proc :',&
392                                     mpi_rank_local,' WARNING lon alloc'
393
394    lnx = size(lat,dim=1)
395    lny = size(lat,dim=2)
396
397    allocate(lat8(lnx,lny),stat=ierror)
398    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
399                                     mpi_rank_local,' WARNING lat alloc'
400
401    lon8 = lon
402    lat8 = lat
403    call oasis_write_grid_r8(cgrid,nx,ny,lon8,lat8,partid=lpartid)
404    deallocate(lon8)
405    deallocate(lat8)
406
407    if (local_timers_on) call oasis_timer_stop(trim(subname))
408    call oasis_debug_exit(subname)
409
410  END SUBROUTINE oasis_write_grid_r4
411
412!--------------------------------------------------------------------------
413
414!> User interface to set angle for 8 byte reals
415
416    SUBROUTINE oasis_write_angle_r8(cgrid, nx, ny, angle, partid)
417
418    !-------------------------------------------------
419    ! Routine to add angles to an existing grid file.
420    !-------------------------------------------------
421
422    implicit none
423
424    character(len=*),         intent (in) :: cgrid       !< grid name
425    integer(kind=ip_intwp_p), intent (in) :: nx          !< global nx size
426    integer(kind=ip_intwp_p), intent (in) :: ny          !< global ny size
427    real(kind=ip_double_p),   intent (in) :: angle(:,:)  !< angles
428    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
429    !-------------------------------------------------
430    integer(kind=ip_intwp_p) :: GRIDID
431    integer(kind=ip_intwp_p) :: ierror
432    integer(kind=ip_intwp_p) :: lnx,lny
433    character(len=*),parameter :: subname = '(oasis_write_angle_r8)'
434    !-------------------------------------------------
435
436    call oasis_debug_enter(subname)
437    if (local_timers_on) call oasis_timer_start(trim(subname))
438
439    if (OASIS_debug >= 15) then
440       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
441    endif
442
443    call oasis_findgrid(cgrid,nx,ny,gridID)
444
445    lnx = size(angle,dim=1)
446    lny = size(angle,dim=2)
447
448    allocate(prism_grid(gridID)%angle(lnx,lny),stat=ierror)
449    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
450                                     mpi_rank_local,' WARNING angle alloc'
451
452    prism_grid(gridID)%angle = angle
453    prism_grid(gridID)%angle_set = .true.
454    if (present(partid)) then
455       if (prism_grid(gridID)%partid > 0 .and. prism_grid(gridID)%partid /= partid) then
456          write(nulprt,*) subname,estr,'partid inconsistency',gridID,prism_grid(gridID)%partid,partid
457          call oasis_abort(file=__FILE__,line=__LINE__)
458       endif
459       prism_grid(gridID)%partid = partid
460       if (OASIS_debug >= 15) then
461          write(nulprt,*) subname,' partid = ',trim(cgrid),partid
462       endif
463    endif
464
465    if (local_timers_on) call oasis_timer_stop(trim(subname))
466    call oasis_debug_exit(subname)
467
468  END SUBROUTINE oasis_write_angle_r8
469
470!--------------------------------------------------------------------------
471
472!> User interface to set angle for 4 byte reals
473
474    SUBROUTINE oasis_write_angle_r4(cgrid, nx, ny, angle, partid)
475
476    !-------------------------------------------------
477    ! Routine to add angles to an existing grid file.
478    !-------------------------------------------------
479
480    implicit none
481
482    character(len=*),         intent (in) :: cgrid       !< grid name
483    integer(kind=ip_intwp_p), intent (in) :: nx          !< global nx size
484    integer(kind=ip_intwp_p), intent (in) :: ny          !< global ny size
485    real(kind=ip_single_p),   intent (in) :: angle(:,:)  !< angles
486    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
487    !-------------------------------------------------
488    real(kind=ip_double_p),allocatable :: angle8(:,:)
489    integer(kind=ip_intwp_p) :: ierror
490    integer(kind=ip_intwp_p) :: lpartid
491    integer(kind=ip_intwp_p) :: lnx,lny
492    character(len=*),parameter :: subname = '(oasis_write_angle_r4)'
493    !-------------------------------------------------
494
495    call oasis_debug_enter(subname)
496    if (local_timers_on) call oasis_timer_start(trim(subname))
497
498    if (OASIS_debug >= 15) then
499       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
500    endif
501
502    lpartid = -1
503    if (present(partid)) then
504       lpartid = partid
505    endif
506    if (OASIS_debug >= 15) then
507       write(nulprt,*) subname,' partid = ',trim(cgrid),lpartid
508    endif
509
510    lnx = size(angle,dim=1)
511    lny = size(angle,dim=2)
512
513    allocate(angle8(lnx,lny),stat=ierror)
514    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
515                                     mpi_rank_local,' WARNING angle8 alloc'
516
517    angle8 = angle
518    call oasis_write_angle_r8(cgrid,nx,ny,angle8,partid=lpartid)
519
520    deallocate(angle8)
521
522    if (local_timers_on) call oasis_timer_stop(trim(subname))
523    call oasis_debug_exit(subname)
524
525  END SUBROUTINE oasis_write_angle_r4
526
527!--------------------------------------------------------------------------
528
529!> User interface to set corner latitudes and longitudes for 8 byte reals
530
531    SUBROUTINE oasis_write_corner_r8(cgrid, nx, ny, nc, clon, clat, partid)
532
533    !-------------------------------------------------
534    ! Routine to add longitudes and latitudes of grid cell corners to an
535    ! existing grids file.
536    !-------------------------------------------------
537
538    implicit none
539
540    character(len=*),         intent (in) :: cgrid  !< grid name
541    integer(kind=ip_intwp_p), intent (in) :: nx     !< global nx size
542    integer(kind=ip_intwp_p), intent (in) :: ny     !< global ny size
543    integer(kind=ip_intwp_p), intent (in) :: nc     !< number of corners per cell
544    real(kind=ip_double_p),   intent (in) :: clon(:,:,:) !< corner longitudes
545    real(kind=ip_double_p),   intent (in) :: clat(:,:,:) !< corner latitudes
546    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
547    !-------------------------------------------------
548    integer(kind=ip_intwp_p) :: GRIDID
549    integer(kind=ip_intwp_p) :: ierror
550    integer(kind=ip_intwp_p) :: lnx,lny
551    character(len=*),parameter :: subname = '(oasis_write_corner_r8)'
552    !-------------------------------------------------
553
554    call oasis_debug_enter(subname)
555    if (local_timers_on) call oasis_timer_start(trim(subname))
556
557    if (OASIS_debug >= 15) then
558       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
559    endif
560
561    call oasis_findgrid(cgrid,nx,ny,gridID)
562
563    lnx = size(clon,dim=1)
564    lny = size(clon,dim=2)
565
566    allocate(prism_grid(gridID)%clon(lnx,lny,nc),stat=ierror)
567    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
568                                     mpi_rank_local,' WARNING clon alloc'
569
570    lnx = size(clat,dim=1)
571    lny = size(clat,dim=2)
572
573    allocate(prism_grid(gridID)%clat(lnx,lny,nc),stat=ierror)
574    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
575                                     mpi_rank_local,' WARNING clat alloc'
576
577    prism_grid(gridID)%nc = nc
578    prism_grid(gridID)%clon = clon
579    prism_grid(gridID)%clat = clat
580    prism_grid(gridID)%corner_set = .true.
581    if (present(partid)) then
582       if (prism_grid(gridID)%partid > 0 .and. prism_grid(gridID)%partid /= partid) then
583          write(nulprt,*) subname,estr,'partid inconsistency',gridID,prism_grid(gridID)%partid,partid
584          call oasis_abort(file=__FILE__,line=__LINE__)
585       endif
586       prism_grid(gridID)%partid = partid
587       if (OASIS_debug >= 15) then
588          write(nulprt,*) subname,' partid = ',trim(cgrid),partid
589       endif
590    endif
591
592    if (local_timers_on) call oasis_timer_stop(trim(subname))
593    call oasis_debug_exit(subname)
594
595  END SUBROUTINE oasis_write_corner_r8
596
597!--------------------------------------------------------------------------
598
599!> User interface to set corner latitudes and longitudes for 4 byte reals
600
601    SUBROUTINE oasis_write_corner_r4(cgrid, nx, ny, nc, clon, clat, partid)
602
603    !-------------------------------------------------
604    ! Routine to add longitudes and latitudes of grid cell corners to an
605    ! existing grids file.
606    !-------------------------------------------------
607
608    implicit none
609
610    character(len=*),         intent (in) :: cgrid  !< grid name
611    integer(kind=ip_intwp_p), intent (in) :: nx     !< global nx size
612    integer(kind=ip_intwp_p), intent (in) :: ny     !< global ny size
613    integer(kind=ip_intwp_p), intent (in) :: nc     !< number of corners per cell
614    real(kind=ip_single_p),   intent (in) :: clon(:,:,:) !< corner longitudes
615    real(kind=ip_single_p),   intent (in) :: clat(:,:,:) !< corner latitudes
616    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
617    !-------------------------------------------------
618    real(kind=ip_double_p), allocatable :: clon8(:,:,:),clat8(:,:,:)
619    integer(kind=ip_intwp_p) :: ierror
620    integer(kind=ip_intwp_p) :: lpartid
621    integer(kind=ip_intwp_p) :: lnx,lny
622    character(len=*),parameter :: subname = '(oasis_write_corner_r4)'
623    !-------------------------------------------------
624
625    call oasis_debug_enter(subname)
626    if (local_timers_on) call oasis_timer_start(trim(subname))
627
628    if (OASIS_debug >= 15) then
629       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
630    endif
631
632    lpartid = -1
633    if (present(partid)) then
634       lpartid = partid
635    endif
636    if (OASIS_debug >= 15) then
637       write(nulprt,*) subname,' partid = ',trim(cgrid),lpartid
638    endif
639
640    lnx = size(clon,dim=1)
641    lny = size(clon,dim=2)
642
643    allocate(clon8(lnx,lny,nc),stat=ierror)
644    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
645                                     mpi_rank_local,' WARNING clon8 alloc'
646
647    lnx = size(clat,dim=1)
648    lny = size(clat,dim=2)
649
650    allocate(clat8(lnx,lny,nc),stat=ierror)
651    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
652                                     mpi_rank_local,' WARNING clat8 alloc'
653
654    clon8 = clon
655    clat8 = clat
656    call oasis_write_corner_r8(cgrid,nx,ny,nc,clon8,clat8,partid=lpartid)
657
658    deallocate(clon8)
659    deallocate(clat8)
660
661    if (local_timers_on) call oasis_timer_stop(trim(subname))
662    call oasis_debug_exit(subname)
663
664  END SUBROUTINE oasis_write_corner_r4
665
666!--------------------------------------------------------------------------
667
668!> User interface to set integer mask values
669
670    SUBROUTINE oasis_write_mask(cgrid, nx, ny, mask, partid, companion)
671
672    !-------------------------------------------------
673    ! Routine to create a new masks file or to add a land see mask to an
674    ! existing masks file.
675    !-------------------------------------------------
676
677    implicit none
678
679    character(len=*),         intent (in) :: cgrid       !< grid name
680    integer(kind=ip_intwp_p), intent (in) :: nx          !< global nx size
681    integer(kind=ip_intwp_p), intent (in) :: ny          !< global ny size
682    integer(kind=ip_intwp_p), intent (in) :: mask(:,:)   !< mask array
683    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
684    character(len=*)        , intent (in),optional :: companion !< companion grid name
685    !-------------------------------------------------
686    integer(kind=ip_intwp_p) :: GRIDID
687    integer(kind=ip_intwp_p) :: ierror
688    integer(kind=ip_intwp_p) :: lnx,lny
689    character(len=*),parameter :: subname = '(oasis_write_mask)'
690    !-------------------------------------------------
691
692    call oasis_debug_enter(subname)
693    if (local_timers_on) call oasis_timer_start(trim(subname))
694
695    if (OASIS_debug >= 15) then
696       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
697    endif
698
699    call oasis_findgrid(cgrid,nx,ny,gridID)
700
701    lnx = size(mask,dim=1)
702    lny = size(mask,dim=2)
703
704    allocate(prism_grid(gridID)%mask(lnx,lny),stat=ierror)
705    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
706                                     mpi_rank_local,' WARNING mask alloc'
707
708    prism_grid(gridID)%mask = mask
709    prism_grid(gridID)%mask_set = .true.
710    if (present(companion)) prism_grid(gridID)%mask_companion = trim(companion)
711    if (present(partid)) then
712       if (prism_grid(gridID)%partid > 0 .and. prism_grid(gridID)%partid /= partid) then
713          write(nulprt,*) subname,estr,'partid inconsistency',gridID,prism_grid(gridID)%partid,partid
714          call oasis_abort(file=__FILE__,line=__LINE__)
715       endif
716       prism_grid(gridID)%partid = partid
717       if (OASIS_debug >= 15) then
718          write(nulprt,*) subname,' partid = ',trim(cgrid),partid
719       endif
720    endif
721
722    if (local_timers_on) call oasis_timer_stop(trim(subname))
723    call oasis_debug_exit(subname)
724
725  END SUBROUTINE oasis_write_mask
726
727!--------------------------------------------------------------------------
728
729!> User interface to set area values for 8 byte reals
730
731    SUBROUTINE oasis_write_area_r8(cgrid, nx, ny, area, partid)
732
733    !-------------------------------------------------
734    ! Routine to create a new areas file or to add areas of a grid to an
735    ! existing areas file.
736    !-------------------------------------------------
737
738    implicit none
739
740    character(len=*),         intent (in) :: cgrid       !< grid name
741    integer(kind=ip_intwp_p), intent (in) :: nx          !< global nx size
742    integer(kind=ip_intwp_p), intent (in) :: ny          !< global ny size
743    real(kind=ip_double_p),   intent (in) :: area(:,:)   !< areas
744    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
745    !-------------------------------------------------
746    integer(kind=ip_intwp_p) :: GRIDID
747    integer(kind=ip_intwp_p) :: ierror
748    integer(kind=ip_intwp_p) :: lnx,lny
749    character(len=*),parameter :: subname = '(oasis_write_area_r8)'
750    !-------------------------------------------------
751
752    call oasis_debug_enter(subname)
753    if (local_timers_on) call oasis_timer_start(trim(subname))
754
755    if (OASIS_debug >= 15) then
756       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
757    endif
758
759    call oasis_findgrid(cgrid,nx,ny,gridID)
760
761    lnx = size(area,dim=1)
762    lny = size(area,dim=2)
763
764    allocate(prism_grid(gridID)%area(lnx,lny),stat=ierror)
765    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
766                                     mpi_rank_local,' WARNING area alloc'
767
768    prism_grid(gridID)%area = area
769    prism_grid(gridID)%area_set = .true.
770    if (present(partid)) then
771       if (prism_grid(gridID)%partid > 0 .and. prism_grid(gridID)%partid /= partid) then
772          write(nulprt,*) subname,estr,'partid inconsistency',gridID,prism_grid(gridID)%partid,partid
773          call oasis_abort(file=__FILE__,line=__LINE__)
774       endif
775       prism_grid(gridID)%partid = partid
776       if (OASIS_debug >= 15) then
777          write(nulprt,*) subname,' partid = ',trim(cgrid),partid
778       endif
779    endif
780
781    if (local_timers_on) call oasis_timer_stop(trim(subname))
782    call oasis_debug_exit(subname)
783
784  END SUBROUTINE oasis_write_area_r8
785
786!--------------------------------------------------------------------------
787
788!> User interface to set area values for 4 byte reals
789
790    SUBROUTINE oasis_write_area_r4(cgrid, nx, ny, area, partid)
791
792    !-------------------------------------------------
793    ! Routine to create a new areas file or to add areas of a grid to an
794    ! existing areas file.
795    !-------------------------------------------------
796
797    implicit none
798
799    character(len=*),         intent (in) :: cgrid       !< grid name
800    integer(kind=ip_intwp_p), intent (in) :: nx          !< global nx size
801    integer(kind=ip_intwp_p), intent (in) :: ny          !< global ny size
802    real(kind=ip_single_p),   intent (in) :: area(:,:)   !< areas
803    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
804    !-------------------------------------------------
805    real(kind=ip_double_p), allocatable :: area8(:,:)
806    integer(kind=ip_intwp_p) :: ierror
807    integer(kind=ip_intwp_p) :: lpartid
808    integer(kind=ip_intwp_p) :: lnx,lny
809    character(len=*),parameter :: subname = '(oasis_write_area_r4)'
810    !-------------------------------------------------
811
812    call oasis_debug_enter(subname)
813    if (local_timers_on) call oasis_timer_start(trim(subname))
814
815    if (OASIS_debug >= 15) then
816       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
817    endif
818
819    lpartid = -1
820    if (present(partid)) then
821       lpartid = partid
822    endif
823    if (OASIS_debug >= 15) then
824       write(nulprt,*) subname,' partid = ',trim(cgrid),lpartid
825    endif
826
827    lnx = size(area,dim=1)
828    lny = size(area,dim=2)
829
830    allocate(area8(lnx,lny),stat=ierror)
831    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
832                                     mpi_rank_local,' WARNING area8 alloc'
833
834    area8 = area
835    call oasis_write_area_r8(cgrid,nx,ny,area8,partid=lpartid)
836
837    deallocate(area8)
838
839    if (local_timers_on) call oasis_timer_stop(trim(subname))
840    call oasis_debug_exit(subname)
841
842  END SUBROUTINE oasis_write_area_r4
843
844!--------------------------------------------------------------------------
845!> User interface to set frac values for 8 byte reals
846
847    SUBROUTINE oasis_write_frac_r8(cgrid, nx, ny, frac, partid, companion)
848
849    !-------------------------------------------------
850    ! Routine to create a new fracs file or to add fracs of a grid to an
851    ! existing fracs file.
852    !-------------------------------------------------
853
854    implicit none
855
856    character(len=*),         intent (in) :: cgrid       !< grid name
857    integer(kind=ip_intwp_p), intent (in) :: nx          !< global nx size
858    integer(kind=ip_intwp_p), intent (in) :: ny          !< global ny size
859    real(kind=ip_double_p),   intent (in) :: frac(:,:)   !< fracs
860    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
861    character(len=*)        , intent (in),optional :: companion !< companion grid name
862    !-------------------------------------------------
863    integer(kind=ip_intwp_p) :: GRIDID
864    integer(kind=ip_intwp_p) :: ierror
865    integer(kind=ip_intwp_p) :: lnx,lny
866    character(len=*),parameter :: subname = '(oasis_write_frac_r8)'
867    !-------------------------------------------------
868
869    call oasis_debug_enter(subname)
870    if (local_timers_on) call oasis_timer_start(trim(subname))
871
872    if (OASIS_debug >= 15) then
873       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
874    endif
875
876    call oasis_findgrid(cgrid,nx,ny,gridID)
877
878    lnx = size(frac,dim=1)
879    lny = size(frac,dim=2)
880
881    allocate(prism_grid(gridID)%frac(lnx,lny),stat=ierror)
882    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
883                                     mpi_rank_local,' WARNING frac alloc'
884
885    prism_grid(gridID)%frac = frac
886    prism_grid(gridID)%frac_set = .true.
887    if (present(companion)) prism_grid(gridID)%frac_companion = trim(companion)
888    if (present(partid)) then
889       if (prism_grid(gridID)%partid > 0 .and. prism_grid(gridID)%partid /= partid) then
890          write(nulprt,*) subname,estr,'partid inconsistency',gridID,prism_grid(gridID)%partid,partid
891          call oasis_abort(file=__FILE__,line=__LINE__)
892       endif
893       prism_grid(gridID)%partid = partid
894       if (OASIS_debug >= 15) then
895          write(nulprt,*) subname,' partid = ',trim(cgrid),partid
896       endif
897    endif
898
899    if (local_timers_on) call oasis_timer_stop(trim(subname))
900    call oasis_debug_exit(subname)
901
902  END SUBROUTINE oasis_write_frac_r8
903
904!--------------------------------------------------------------------------
905
906!> User interface to set frac values for 4 byte reals
907
908    SUBROUTINE oasis_write_frac_r4(cgrid, nx, ny, frac, partid, companion)
909
910    !-------------------------------------------------
911    ! Routine to create a new fracs file or to add fracs of a grid to an
912    ! existing fracs file.
913    !-------------------------------------------------
914
915    implicit none
916
917    character(len=*),         intent (in) :: cgrid       !< grid name
918    integer(kind=ip_intwp_p), intent (in) :: nx          !< global nx size
919    integer(kind=ip_intwp_p), intent (in) :: ny          !< global ny size
920    real(kind=ip_single_p),   intent (in) :: frac(:,:)   !< fracs
921    integer(kind=ip_intwp_p), intent (in),optional :: partid  !< partition id if nonglobal data
922    character(len=*)        , intent (in),optional :: companion !< companion grid name
923    !-------------------------------------------------
924    real(kind=ip_double_p), allocatable :: frac8(:,:)
925    integer(kind=ip_intwp_p) :: ierror
926    integer(kind=ip_intwp_p) :: lpartid
927    integer(kind=ip_intwp_p) :: lnx,lny
928    character(len=*),parameter :: subname = '(oasis_write_frac_r4)'
929    !-------------------------------------------------
930
931    call oasis_debug_enter(subname)
932    if (local_timers_on) call oasis_timer_start(trim(subname))
933
934    if (OASIS_debug >= 15) then
935       write(nulprt,*) subname,' size = ',trim(cgrid),nx,ny
936    endif
937
938    lpartid = -1
939    if (present(partid)) then
940       lpartid = partid
941    endif
942    if (OASIS_debug >= 15) then
943       write(nulprt,*) subname,' partid = ',trim(cgrid),lpartid
944    endif
945
946    lnx = size(frac,dim=1)
947    lny = size(frac,dim=2)
948
949    allocate(frac8(lnx,lny),stat=ierror)
950    if (ierror /= 0) write(nulprt,*) subname,' model :',compid,' proc :',&
951                                     mpi_rank_local,' WARNING frac8 alloc'
952
953    frac8 = frac
954    if (present(companion)) then
955       call oasis_write_frac_r8(cgrid,nx,ny,frac8,partid=lpartid,companion=trim(companion))
956    else
957       call oasis_write_frac_r8(cgrid,nx,ny,frac8,partid=lpartid)
958    endif
959
960    deallocate(frac8)
961
962    if (local_timers_on) call oasis_timer_stop(trim(subname))
963    call oasis_debug_exit(subname)
964
965  END SUBROUTINE oasis_write_frac_r4
966
967!--------------------------------------------------------------------------
968
969!> User interface to indicate user defined grids are done
970
971    SUBROUTINE oasis_terminate_grids_writing()
972    !-------------------------------------------------
973    ! Routine to terminate the grids writing.
974    !-------------------------------------------------
975
976    implicit none
977    integer(kind=ip_i4_p) :: n
978    character(len=*),parameter :: subname = '(oasis_terminate_grids_writing)'
979
980    call oasis_debug_enter(subname)
981    if (local_timers_on) call oasis_timer_start(trim(subname))
982
983    if (OASIS_debug >= 15) then
984       write(nulprt,*) subname,' prism_ngrid = ',prism_ngrid
985    endif
986
987    do n = 1,prism_ngrid
988       prism_grid(n)%terminated = .true.
989    enddo
990
991    call oasis_print_grid_data()
992! moved to prism_method_enddef for synchronization
993!    call oasis_write2files()
994
995    if (local_timers_on) call oasis_timer_stop(trim(subname))
996    call oasis_debug_exit(subname)
997
998  END SUBROUTINE oasis_terminate_grids_writing
999
1000!--------------------------------------------------------------------------
1001
1002!> Interface that actually writes fields to grid files
1003
1004    SUBROUTINE oasis_write2files()
1005
1006    !-------------------------------------------------
1007    !> Write fields to grid files.
1008    !> Only write fields that have been buffered and
1009    !> if prism_grid_terminate_grids_writing has been called.
1010    !> This is called by all tasks from oasis_enddef.
1011    !-------------------------------------------------
1012
1013    implicit none
1014
1015    !-------------------------------------------------
1016    character(len=ic_med) :: filename  ! grid filename
1017    character(len=ic_med) :: fldname   ! full field name
1018    character(len=ic_med) :: cgrid     ! grid name
1019    logical :: exists                  ! check if file exists
1020    integer(kind=ip_i4_p) :: m,n,n1,g,p  ! counter
1021    integer(kind=ip_i4_p) :: partid    ! part id
1022    integer(kind=ip_i4_p) :: taskid    ! task id for writing
1023    integer(kind=ip_i4_p) :: nx,ny,nc  ! grid size
1024    integer(kind=ip_i4_p) :: tnx,tny   ! temporary grid size
1025    logical :: partid_grid             ! global or local grid
1026    logical :: active_task             ! determine which tasks are involved
1027    logical :: write_task              ! task for writing
1028    real(kind=ip_realwp_p),allocatable :: rloc(:,:) ! local array
1029    real(kind=ip_realwp_p),allocatable :: rglo(:,:) ! global array
1030    real(kind=ip_realwp_p),allocatable :: r3glo(:,:,:) ! global array
1031    integer(kind=ip_i4_p) ,allocatable :: iglo(:,:) ! global array
1032    integer(kind=ip_intwp_p) :: gcnt
1033    logical                  :: found
1034    character(len=ic_med)   ,allocatable :: gname0(:),gname(:)
1035    character(len=ic_lvar2) ,allocatable :: pname0(:),pname(:)
1036    character(len=*),parameter :: undefined_partname = '(UnDeFiNeD_PArtnaME)'
1037    character(len=*),parameter :: subname = '(oasis_write2files)'
1038    !-------------------------------------------------
1039
1040    call oasis_debug_enter(subname)
1041    call oasis_timer_start('oasis_gridw')
1042
1043    call oasis_mpi_bcast(writing_grids_call,mpi_comm_local,subname//'writing_grids_call')
1044    if (writing_grids_call .eq. 1) then
1045
1046    if (local_timers_on) call oasis_timer_start('oasis_gridw_reducelists')
1047    allocate(gname0(prism_ngrid))
1048    allocate(pname0(prism_ngrid))
1049    do n = 1,prism_ngrid
1050       gname0(n) = prism_grid(n)%gridname
1051       if (prism_grid(n)%partid > 0 .and. prism_grid(n)%partid <= prism_npart) then
1052          pname0(n) = prism_part(prism_grid(n)%partid)%partname
1053       elseif (prism_grid(n)%partid == -1) then
1054          pname0(n) = undefined_partname
1055       else
1056          write(nulprt,*) subname,estr,'illegal partition id for grid ',trim(prism_grid(n)%gridname),prism_grid(n)%partid
1057       endif
1058    enddo
1059
1060    call oasis_mpi_reducelists(gname0,mpi_comm_local,gcnt,gname,'write2files',fastcheck=.true., &
1061         linp2=pname0,lout2=pname,spval2=undefined_partname)
1062    deallocate(gname0)
1063    deallocate(pname0)
1064    if (local_timers_on) call oasis_timer_stop('oasis_gridw_reducelists')
1065
1066    !-------------------------------------
1067    !> * Check that a grid defined on a partitition is defined on all tasks on that partition.
1068    !-------------------------------------
1069
1070    if (local_timers_on) call oasis_timer_start('oasis_gridw_chkpart')
1071    do n = 1,gcnt
1072       if (pname(n) /= undefined_partname) then
1073          do p = 1,prism_npart
1074             if (pname(n) == prism_part(p)%partname .and. prism_part(p)%mpicom /= MPI_COMM_NULL) then
1075                found = .false.
1076                do g = 1,prism_ngrid
1077                   if (prism_grid(g)%gridname == gname(n)) found = .true.
1078                enddo
1079                if (.not. found) then
1080                   write(nulprt,*) subname,estr,'grid with partition not defined on all partition tasks: ',trim(gname(n))
1081                   call oasis_abort(file=__FILE__,line=__LINE__)
1082                endif
1083             endif
1084          enddo
1085       endif
1086    enddo
1087    if (local_timers_on) call oasis_timer_stop('oasis_gridw_chkpart')
1088
1089    if (local_timers_on) call oasis_timer_start('oasis_gridw_writefiles')
1090
1091    !-------------------------------------
1092    !> * Write grid information
1093    !-------------------------------------
1094
1095    do g = 1,gcnt
1096    do n = 1,prism_ngrid
1097    if (prism_grid(n)%terminated) then
1098    if (prism_grid(n)%gridname == gname(g)) then
1099
1100       cgrid  = gname(g)
1101       partid = prism_grid(n)%partid
1102       prism_grid(n)%written = .true.
1103       tnx = -1
1104       tny = -1
1105
1106       !-------------------------------------
1107       !>   * Determine which tasks are associated with the grid information
1108       !-------------------------------------
1109
1110       active_task = .false.
1111       write_task = .false.
1112       if (pname(g) == undefined_partname) then
1113          partid_grid = .false.
1114          if (mpi_rank_local == 0) active_task = .true.
1115          if (mpi_rank_local == 0) write_task = .true.
1116       else
1117          partid_grid = .true.
1118          if (partid > 0 .and. partid <= prism_npart) then
1119             taskid = 0
1120             if (prism_part(partid)%mpicom /= MPI_COMM_NULL) active_task = .true.
1121             if (prism_part(partid)%rank == taskid) write_task = .true.
1122          elseif (partid == -1) then
1123             active_task = .false.
1124             write_task = .false.
1125          else
1126             write(nulprt,*) subname,estr,'illegal partid for grid:',trim(gname(g)),trim(pname(g)),partid
1127             call oasis_abort(file=__FILE__,line=__LINE__)
1128          endif
1129       endif
1130
1131       if (OASIS_debug >= 15) then
1132          write(nulprt,*) subname,' ',trim(gname(g)),':',trim(pname(g)),': partid_grid=', &
1133                          partid_grid,'active_task=',active_task,'write_task=',write_task
1134       endif
1135
1136       if (active_task) then
1137
1138         nx = prism_grid(n)%nx
1139         ny = prism_grid(n)%ny
1140         nc = prism_grid(n)%nc
1141
1142         allocate(rglo(nx,ny))
1143
1144         !-------------------------------------
1145         !>   * Check that array sizes match for all fields
1146         !-------------------------------------
1147
1148         if (prism_grid(n)%grid_set) then
1149           if (tnx <= 0 .or. tny <= 0) then
1150             tnx = size(prism_grid(n)%lon,dim=1)
1151             tny = size(prism_grid(n)%lon,dim=2)
1152           endif
1153           if (size(prism_grid(n)%lon,dim=1) /= tnx .or. &
1154               size(prism_grid(n)%lon,dim=2) /= tny .or. &
1155               size(prism_grid(n)%lat,dim=1) /= tnx .or. &
1156               size(prism_grid(n)%lat,dim=2) /= tny ) then
1157             write(nulprt,*) subname,estr,'inconsistent array size lon/lat ',tnx,tny, &
1158               size(prism_grid(n)%lon,dim=1),size(prism_grid(n)%lon,dim=2),  &
1159               size(prism_grid(n)%lat,dim=1),size(prism_grid(n)%lat,dim=2)
1160             call oasis_abort(file=__FILE__,line=__LINE__)
1161           endif
1162
1163           !-------------------------------------
1164           !>   * Gather longitudes if needed and write from root
1165           !-------------------------------------
1166
1167           filename = 'grids.nc'
1168           fldname  = trim(cgrid)//'.lon'
1169           if (.not.oasis_io_varexists(filename,fldname)) then
1170             if (partid_grid) then
1171               call oasis_grid_loc2glo(prism_grid(n)%lon,rglo,partid,0)
1172             else
1173               rglo = prism_grid(n)%lon
1174             endif
1175             if (write_task) call oasis_io_write_2dgridfld_fromroot(filename,fldname,rglo,nx,ny)
1176           endif
1177
1178           !-------------------------------------
1179           !>   * Gather latitudes if needed and write from root
1180           !-------------------------------------
1181
1182           filename = 'grids.nc'
1183           fldname  = trim(cgrid)//'.lat'
1184           if (.not.oasis_io_varexists(filename,fldname)) then
1185             if (partid_grid) then
1186               call oasis_grid_loc2glo(prism_grid(n)%lat,rglo,partid,0)
1187             else
1188               rglo = prism_grid(n)%lat
1189             endif
1190             if (write_task) call oasis_io_write_2dgridfld_fromroot(filename,fldname,rglo,nx,ny)
1191           endif
1192         endif  ! grid_set
1193
1194         if (prism_grid(n)%corner_set) then
1195           if (tnx <= 0 .or. tny <= 0) then
1196             tnx = size(prism_grid(n)%clon,dim=1)
1197             tny = size(prism_grid(n)%clon,dim=2)
1198           endif
1199           if (size(prism_grid(n)%clon,dim=1) /= tnx .or. &
1200               size(prism_grid(n)%clon,dim=2) /= tny .or. &
1201               size(prism_grid(n)%clat,dim=1) /= tnx .or. &
1202               size(prism_grid(n)%clat,dim=2) /= tny ) then
1203             write(nulprt,*) subname,estr,'inconsistent array size clon/clat ',tnx,tny, &
1204               size(prism_grid(n)%clon,dim=1),size(prism_grid(n)%clon,dim=2),  &
1205               size(prism_grid(n)%clat,dim=1),size(prism_grid(n)%clat,dim=2)
1206             call oasis_abort(file=__FILE__,line=__LINE__)
1207           endif
1208
1209           !-------------------------------------
1210           !>   * Gather corner longitudes if needed and write from root
1211           !-------------------------------------
1212
1213           allocate(r3glo(nx,ny,nc))
1214           filename = 'grids.nc'
1215           fldname  = trim(cgrid)//'.clo'
1216           if (.not.oasis_io_varexists(filename,fldname)) then
1217             if (partid_grid) then
1218               allocate(rloc(tnx,tny))
1219               do n1 = 1,nc
1220                 rloc(:,:) = prism_grid(n)%clon(:,:,n1)
1221                 call oasis_grid_loc2glo(rloc,rglo,partid,0)
1222                 r3glo(:,:,n1) = rglo(:,:)
1223               enddo
1224               deallocate(rloc)
1225             else
1226               r3glo = prism_grid(n)%clon
1227             endif
1228             if (write_task) call oasis_io_write_3dgridfld_fromroot(filename,fldname,r3glo,nx,ny,nc)
1229           endif
1230
1231           !-------------------------------------
1232           !>   * Gather corner latitudes if needed and write from root
1233           !-------------------------------------
1234
1235           filename = 'grids.nc'
1236           fldname  = trim(cgrid)//'.cla'
1237           if (.not.oasis_io_varexists(filename,fldname)) then
1238             if (partid_grid) then
1239               allocate(rloc(tnx,tny))
1240               do n1 = 1,nc
1241                 rloc(:,:) = prism_grid(n)%clat(:,:,n1)
1242                 call oasis_grid_loc2glo(rloc,rglo,partid,0)
1243                 r3glo(:,:,n1) = rglo(:,:)
1244               enddo
1245               deallocate(rloc)
1246             else
1247               r3glo = prism_grid(n)%clat
1248             endif
1249             if (write_task) call oasis_io_write_3dgridfld_fromroot(filename,fldname,r3glo,nx,ny,nc)
1250           endif
1251           deallocate(r3glo)
1252         endif  ! corner_set
1253
1254         if (prism_grid(n)%area_set) then
1255           if (tnx <= 0 .or. tny <= 0) then
1256             tnx = size(prism_grid(n)%area,dim=1)
1257             tny = size(prism_grid(n)%area,dim=2)
1258           endif
1259           if (size(prism_grid(n)%area,dim=1) /= tnx .or. &
1260               size(prism_grid(n)%area,dim=2) /= tny ) then
1261             write(nulprt,*) subname,estr,'inconsistent array size area ',tnx,tny, &
1262               size(prism_grid(n)%area,dim=1),size(prism_grid(n)%area,dim=2)
1263             call oasis_abort(file=__FILE__,line=__LINE__)
1264           endif
1265
1266           !-------------------------------------
1267           !>   * Gather areas if needed and write from root
1268           !-------------------------------------
1269
1270           filename = 'areas.nc'
1271           fldname  = trim(cgrid)//'.srf'
1272           if (.not.oasis_io_varexists(filename,fldname)) then
1273             if (partid_grid) then
1274               call oasis_grid_loc2glo(prism_grid(n)%area,rglo,partid,0)
1275             else
1276               rglo = prism_grid(n)%area
1277             endif
1278             if (write_task) call oasis_io_write_2dgridfld_fromroot(filename,fldname,rglo,nx,ny)
1279           endif
1280         endif  ! area_set
1281
1282         if (prism_grid(n)%frac_set) then
1283           if (tnx <= 0 .or. tny <= 0) then
1284             tnx = size(prism_grid(n)%frac,dim=1)
1285             tny = size(prism_grid(n)%frac,dim=2)
1286           endif
1287           if (size(prism_grid(n)%frac,dim=1) /= tnx .or. &
1288               size(prism_grid(n)%frac,dim=2) /= tny ) then
1289             write(nulprt,*) subname,estr,'inconsistent array size frac ',tnx,tny, &
1290               size(prism_grid(n)%frac,dim=1),size(prism_grid(n)%frac,dim=2)
1291             call oasis_abort(file=__FILE__,line=__LINE__)
1292           endif
1293
1294           !-------------------------------------
1295           !>   * Gather fracs if needed and write from root
1296           !-------------------------------------
1297
1298           filename = 'masks.nc'
1299           fldname  = trim(cgrid)//'.frc'
1300           if (.not.oasis_io_varexists(filename,fldname)) then
1301             if (partid_grid) then
1302               call oasis_grid_loc2glo(prism_grid(n)%frac,rglo,partid,0)
1303             else
1304               rglo = prism_grid(n)%frac
1305             endif
1306             if (write_task) then
1307                call oasis_io_write_2dgridfld_fromroot(filename,fldname,rglo,nx,ny)
1308                call oasis_io_write_fldattr(filename,fldname,'coherent_with_grid',trim(prism_grid(n)%frac_companion))
1309             endif
1310           endif
1311         endif  ! frac_set
1312
1313         if (prism_grid(n)%angle_set) then
1314           if (tnx <= 0 .or. tny <= 0) then
1315             tnx = size(prism_grid(n)%angle,dim=1)
1316             tny = size(prism_grid(n)%angle,dim=2)
1317           endif
1318           if (size(prism_grid(n)%angle,dim=1) /= tnx .or. &
1319               size(prism_grid(n)%angle,dim=2) /= tny ) then
1320             write(nulprt,*) subname,estr,'inconsistent array size angle ',tnx,tny, &
1321                size(prism_grid(n)%angle,dim=1),size(prism_grid(n)%angle,dim=2)
1322             call oasis_abort(file=__FILE__,line=__LINE__)
1323           endif
1324
1325           !-------------------------------------
1326           !>   * Gather angles if needed and write from root
1327           !-------------------------------------
1328
1329           filename = 'grids.nc'
1330           fldname  = trim(cgrid)//'.ang'
1331           if (.not.oasis_io_varexists(filename,fldname)) then
1332             if (partid_grid) then
1333               call oasis_grid_loc2glo(prism_grid(n)%lon,rglo,partid,0)
1334             else
1335               rglo = prism_grid(n)%angle
1336             endif
1337             if (write_task) call oasis_io_write_2dgridfld_fromroot(filename,fldname,rglo,nx,ny)
1338           endif
1339         endif  ! angle_set
1340
1341         if (prism_grid(n)%mask_set) then
1342           if (tnx <= 0 .or. tny <= 0) then
1343             tnx = size(prism_grid(n)%mask,dim=1)
1344             tny = size(prism_grid(n)%mask,dim=2)
1345           endif
1346           if (size(prism_grid(n)%mask,dim=1) /= tnx .or. &
1347               size(prism_grid(n)%mask,dim=2) /= tny ) then
1348             write(nulprt,*) subname,estr,'inconsistent array size mask ',tnx,tny, &
1349                size(prism_grid(n)%mask,dim=1),size(prism_grid(n)%mask,dim=2)
1350             call oasis_abort(file=__FILE__,line=__LINE__)
1351           endif
1352
1353           !-------------------------------------
1354           !>   * Gather masks if needed and write from root
1355           !-------------------------------------
1356
1357           filename = 'masks.nc'
1358           fldname  = trim(cgrid)//'.msk'
1359           if (.not.oasis_io_varexists(filename,fldname)) then
1360             allocate(iglo(nx,ny))
1361             if (partid_grid) then
1362               allocate(rloc(tnx,tny))
1363               rloc = prism_grid(n)%mask
1364               call oasis_grid_loc2glo(rloc,rglo,partid,0)
1365               iglo = nint(rglo)
1366               deallocate(rloc)
1367             else
1368               iglo = prism_grid(n)%mask
1369             endif
1370             if (write_task) then
1371                call oasis_io_write_2dgridint_fromroot(filename,fldname,iglo,nx,ny)
1372                call oasis_io_write_fldattr(filename,fldname,'coherent_with_grid',trim(prism_grid(n)%mask_companion))
1373             endif
1374             deallocate(iglo)
1375           endif
1376         endif ! mask_set
1377
1378         deallocate(rglo)
1379
1380       endif  ! active_task
1381
1382    endif  ! grid_name = gname
1383    endif  ! terminated
1384    enddo  ! n = 1,prism_ngrid
1385    enddo  ! g = 1,gcnt
1386
1387    deallocate(gname,pname)
1388
1389    if (local_timers_on) call oasis_timer_stop('oasis_gridw_writefiles')
1390    endif ! writing_grids_call
1391
1392    call oasis_timer_stop('oasis_gridw')
1393    call oasis_debug_exit(subname)
1394
1395  END SUBROUTINE oasis_write2files
1396!--------------------------------------------------------------------------
1397
1398!> Local interface to find gridID for a specified grid name
1399
1400    SUBROUTINE oasis_findgrid(cgrid,nx,ny,gridID)
1401    !-------------------------------------------------
1402    ! Routine that sets gridID, identifies existing
1403    ! grid with cgrid name or starts a new one
1404    !-------------------------------------------------
1405    implicit none
1406
1407    character(len=*),         intent (in) :: cgrid       !< grid name
1408    integer(kind=ip_intwp_p), intent (in) :: nx          !< global nx size
1409    integer(kind=ip_intwp_p), intent (in) :: ny          !< global ny size
1410    integer(kind=ip_intwp_p), intent(out) :: gridID      !< gridID matching cgrid
1411    !-------------------------------------------------
1412    integer(kind=ip_intwp_p) :: n
1413    character(len=*),parameter :: subname = '(oasis_findgrid)'
1414    !-------------------------------------------------
1415
1416    call oasis_debug_enter(subname)
1417    if (local_timers_on) call oasis_timer_start(trim(subname))
1418
1419    gridID = -1
1420    do n = 1,prism_ngrid
1421       if (trim(cgrid) == trim(prism_grid(n)%gridname)) then
1422          gridID = n
1423          ! since grid is defined before, make sure nx,ny match
1424          if (nx /= prism_grid(gridID)%nx .or. ny /= prism_grid(gridID)%ny) then
1425             write(nulprt,*) subname,estr,'in predefined grid size = ',nx,ny, &
1426                prism_grid(gridID)%nx,prism_grid(gridID)%ny
1427             call oasis_abort(file=__FILE__,line=__LINE__)
1428          endif
1429       endif
1430    enddo
1431
1432    if (gridID < 1) then
1433       prism_ngrid = prism_ngrid+1
1434       gridID = prism_ngrid
1435    endif
1436
1437    prism_grid(gridID)%gridname = trim(cgrid)
1438    prism_grid(gridID)%nx = nx
1439    prism_grid(gridID)%ny = ny
1440
1441    if (local_timers_on) call oasis_timer_stop(trim(subname))
1442    call oasis_debug_exit(subname)
1443
1444  END SUBROUTINE oasis_findgrid
1445
1446!--------------------------------------------------------------------------
1447
1448!> Local routine that gathers the local array using partition information
1449
1450  SUBROUTINE oasis_grid_loc2glo(aloc,aglo,partid,taskid)
1451    implicit none
1452
1453    real(kind=ip_realwp_p),intent(in)    :: aloc(:,:)  !< local array
1454    real(kind=ip_realwp_p),intent(inout) :: aglo(:,:)  !< global array
1455    integer(kind=ip_i4_p) ,intent(in)    :: partid     !< partition id for local data
1456    integer(kind=ip_i4_p) ,intent(in)    :: taskid     !< task id to gather data to
1457    !-------------------------------------------------
1458    type(mct_aVect) :: avloc,avglo
1459    integer(kind=ip_i4_p) :: i,j,n
1460    integer(kind=ip_i4_p) :: lnx,lny,gnx,gny
1461    character(len=*),parameter :: subname = '(oasis_grid_loc2glo)'
1462    !-------------------------------------------------
1463
1464    call oasis_debug_enter(subname)
1465    if (local_timers_on) call oasis_timer_start(trim(subname))
1466
1467    if (prism_part(partid)%mpicom /= MPI_COMM_NULL) then
1468
1469       lnx = size(aloc,dim=1)
1470       lny = size(aloc,dim=2)
1471       gnx = size(aglo,dim=1)
1472       gny = size(aglo,dim=2)
1473       call mct_avect_init(avloc,rList='field',lsize=lnx*lny)
1474
1475       n = 0
1476       do j = 1,lny
1477       do i = 1,lnx
1478          n = n + 1
1479          avloc%rattr(1,n) = aloc(i,j)
1480       enddo
1481       enddo
1482
1483       call mct_aVect_gather(avloc,avglo,prism_part(partid)%pgsmap,taskid,prism_part(partid)%mpicom)
1484
1485       if (prism_part(partid)%rank == taskid) then
1486          n = 0
1487          do j = 1,gny
1488          do i = 1,gnx
1489             n = n + 1
1490             aglo(i,j) = avglo%rattr(1,n)
1491          enddo
1492          enddo
1493          call mct_avect_clean(avglo)
1494       endif
1495
1496       call mct_avect_clean(avloc)
1497
1498    endif
1499
1500    if (local_timers_on) call oasis_timer_stop(trim(subname))
1501    call oasis_debug_exit(subname)
1502
1503  END SUBROUTINE oasis_grid_loc2glo
1504
1505!--------------------------------------------------------------------------
1506END MODULE mod_oasis_grid
1507!--------------------------------------------------------------------------
1508
1509
1510     
Note: See TracBrowser for help on using the repository browser.