ocnice_prep  1.13.0
All Data Structures Files Functions Variables Pages
utils_mod.F90
Go to the documentation of this file.
1 
7 module utils_mod
8 
9  use netcdf
10  use arrays_mod, only : eta
11  use init_mod, only : debug, logunit, vardefs, fsrc
12 
13  implicit none
14 
15  private
16 
17  interface getfield
18  module procedure getfield2d
19  module procedure getfield3d
20  end interface getfield
21 
22  interface packarrays
23  module procedure packarrays2d
24  module procedure packarrays3d
25  end interface packarrays
26 
27  interface remap
28  module procedure remap1d
29  module procedure remap2d
30  module procedure remap3d
31  end interface remap
32 
33  interface getvecpair
34  module procedure getvecpair2d
35  module procedure getvecpair3d
36  end interface getvecpair
37 
38  interface dumpnc
39  module procedure dumpnc1d
40  module procedure dumpnc2d
41  module procedure dumpnc3d
42  module procedure dumpnc3dk
43  end interface dumpnc
44 
45  public :: getfield
46  public :: packarrays
47  public :: remap
48  public :: dumpnc
49  public :: nf90_err
50 
51 contains
64  subroutine packarrays2d(filesrc, wgtsdir, cosrot, sinrot, vars, dims, nflds, fields)
65 
66  character(len=*), intent(in) :: filesrc,wgtsdir
67  real(kind=8), intent(in) :: cosrot(:),sinrot(:)
68  type(vardefs), intent(in) :: vars(:)
69  integer, intent(in) :: dims(:)
70  integer, intent(in) :: nflds
71  real(kind=8), intent(out) :: fields(:,:)
72 
73  ! local variables
74  integer :: n, nn
75  real(kind=8), allocatable :: vecpair(:,:)
76  character(len=20) :: subname = 'packarrays2d'
77  !----------------------------------------------------------------------------
78 
79  if (debug)write(logunit,'(a)')'enter '//trim(subname)
80 
81  fields=0.0
82  ! obtain vector pairs
83  do n = 1,nflds
84  if (trim(vars(n)%var_grid) == 'Cu' .or. trim(vars(n)%var_grid) == 'Bu_x') then
85  allocate(vecpair(dims(1)*dims(2),2)); vecpair = 0.0
86  call getvecpair(trim(filesrc), trim(wgtsdir), cosrot, sinrot, &
87  trim(vars(n)%var_name), trim(vars(n)%var_grid(1:2)), &
88  trim(vars(n)%var_pair), trim(vars(n)%var_pair_grid(1:2)), &
89  dims=(/dims(1),dims(2)/), vecpair=vecpair)
90  end if
91  end do
92 
93  ! create packed array
94  nn = 0
95  do n = 1,nflds
96  if (len_trim(vars(n)%var_pair) == 0) then
97  nn = nn + 1
98  call getfield(trim(filesrc), trim(vars(n)%var_name), dims=(/dims(1),dims(2)/), &
99  field=fields(nn,:))
100  else ! fill with vector pairs
101  nn = nn+1
102  ! ocn vectors
103  if (trim(vars(n)%var_grid) == 'Cu')fields(nn,:) = vecpair(:,1)
104  if (trim(vars(n)%var_grid) == 'Cv')fields(nn,:) = vecpair(:,2)
105  ! ice vectors
106  if (trim(vars(n)%var_grid) == 'Bu_x')fields(nn,:) = vecpair(:,1)
107  if (trim(vars(n)%var_grid) == 'Bu_y')fields(nn,:) = vecpair(:,2)
108  end if
109  end do
110 
111  if (debug)write(logunit,'(a)')'exit '//trim(subname)
112  end subroutine packarrays2d
113 
126  subroutine packarrays3d(filesrc, wgtsdir, cosrot, sinrot, vars, dims, nflds, fields)
128  character(len=*), intent(in) :: filesrc,wgtsdir
129  real(kind=8), intent(in) :: cosrot(:),sinrot(:)
130  type(vardefs), intent(in) :: vars(:)
131  integer, intent(in) :: dims(:)
132  integer, intent(in) :: nflds
133  real(kind=8), intent(out) :: fields(:,:,:)
134 
135  ! local variables
136  integer :: n, nn
137  real(kind=8), allocatable :: vecpair(:,:,:)
138  character(len=20) :: subname = 'packarrays3d'
139  !----------------------------------------------------------------------------
140 
141  if (debug)write(logunit,'(a)')'enter '//trim(subname)
142 
143  fields=0.0
144  ! obtain vector pairs
145  do n = 1,nflds
146  if (trim(vars(n)%var_grid) == 'Cu') then
147  allocate(vecpair(dims(3),dims(1)*dims(2),2)); vecpair = 0.0
148  call getvecpair(trim(filesrc), trim(wgtsdir), cosrot, sinrot, &
149  trim(vars(n)%var_name), trim(vars(n)%var_grid), &
150  trim(vars(n)%var_pair), trim(vars(n)%var_pair_grid), &
151  dims=(/dims(1),dims(2),dims(3)/), vecpair=vecpair)
152  end if
153  end do
154 
155  ! create packed array
156  nn = 0
157  do n = 1,nflds
158  if (len_trim(vars(n)%var_pair) == 0) then
159  nn = nn + 1
160  if (trim(vars(n)%var_name) .eq. 'eta') then
161  fields(nn,:,:) = eta(:,:)
162  else
163  call getfield(trim(filesrc), trim(vars(n)%var_name), dims=(/dims(1),dims(2),dims(3)/), &
164  field=fields(nn,:,:))
165  end if
166  else ! fill with vector pairs
167  nn = nn+1
168  if (trim(vars(n)%var_grid) == 'Cu')fields(nn,:,:) = vecpair(:,:,1)
169  if (trim(vars(n)%var_grid) == 'Cv')fields(nn,:,:) = vecpair(:,:,2)
170  end if
171  end do
172 
173  if (debug)write(logunit,'(a)')'exit '//trim(subname)
174  end subroutine packarrays3d
175 
190  subroutine getvecpair2d(fname, wdir, cosrot, sinrot, vname1, vgrid1, vname2, vgrid2, dims, vecpair)
192  character(len=*), intent(in) :: fname
193  character(len=*), intent(in) :: wdir
194  real(kind=8), intent(in) :: cosrot(:), sinrot(:)
195  character(len=*), intent(in) :: vname1, vgrid1, vname2, vgrid2
196  integer, intent(in) :: dims(:)
197  real(kind=8), intent(out) :: vecpair(:,:)
198 
199  ! local variables
200  integer :: ii
201  real(kind=8), dimension(dims(1)*dims(2)) :: urot, vrot
202  character(len=240) :: wgtsfile
203  character(len=20) :: subname = 'getvecpair2d'
204  !----------------------------------------------------------------------------
205 
206  if (debug)write(logunit,'(a)')'enter '//trim(subname)
207 
208  wgtsfile = trim(wdir)//'tripole.'//trim(fsrc)//'.'//vgrid1//'.to.Ct.bilinear.nc'
209  call getfield(fname, vname1, dims=dims, field=vecpair(:,1), wgts=trim(wgtsfile))
210  if (debug)write(logunit,'(a)')'wgtsfile for 2d vector '//trim(vname1)//' '//trim(wgtsfile)
211  wgtsfile = trim(wdir)//'tripole.'//trim(fsrc)//'.'//vgrid2//'.to.Ct.bilinear.nc'
212  call getfield(fname, vname2, dims=dims, field=vecpair(:,2), wgts=trim(wgtsfile))
213  if (debug)write(logunit,'(a)')'wgtsfile for 2d vector '//trim(vname2)//' '//trim(wgtsfile)
214 
215  urot = 0.0; vrot = 0.0
216  do ii = 1,dims(1)*dims(2)
217  urot(ii) = vecpair(ii,1)*cosrot(ii) + vecpair(ii,2)*sinrot(ii)
218  vrot(ii) = vecpair(ii,2)*cosrot(ii) - vecpair(ii,1)*sinrot(ii)
219  end do
220  vecpair(:,1) = urot(:)
221  vecpair(:,2) = vrot(:)
222 
223  if (debug) write(logunit,'(a)')'exit '//trim(subname)
224  end subroutine getvecpair2d
225 
240  subroutine getvecpair3d(fname, wdir, cosrot, sinrot, vname1, vgrid1, vname2, vgrid2, dims, vecpair)
242  character(len=*), intent(in) :: fname
243  character(len=*), intent(in) :: wdir
244  real(kind=8), intent(in) :: cosrot(:), sinrot(:)
245  character(len=*), intent(in) :: vname1, vgrid1, vname2, vgrid2
246  integer, intent(in) :: dims(:)
247  real(kind=8), intent(out) :: vecpair(:,:,:)
248 
249  ! local variables
250  integer :: ii,k
251  real(kind=8), dimension(dims(1)*dims(2)) :: urot, vrot
252  character(len=240) :: wgtsfile
253  character(len=20) :: subname = 'getvecpair3d'
254  !----------------------------------------------------------------------------
255 
256  if (debug)write(logunit,'(a)')'enter '//trim(subname)
257 
258  wgtsfile = trim(wdir)//'tripole.'//trim(fsrc)//'.'//vgrid1//'.to.Ct.bilinear.nc'
259  call getfield(fname, vname1, dims=dims, field=vecpair(:,:,1), wgts=trim(wgtsfile))
260  wgtsfile = trim(wdir)//'tripole.'//trim(fsrc)//'.'//vgrid2//'.to.Ct.bilinear.nc'
261  call getfield(fname, vname2, dims=dims, field=vecpair(:,:,2), wgts=trim(wgtsfile))
262 
263  do k = 1,dims(3)
264  urot = 0.0; vrot = 0.0
265  do ii = 1,dims(1)*dims(2)
266  urot(ii)= vecpair(k,ii,1)*cosrot(ii) + vecpair(k,ii,2)*sinrot(ii)
267  vrot(ii)= vecpair(k,ii,2)*cosrot(ii) - vecpair(k,ii,1)*sinrot(ii)
268  end do
269  vecpair(k,:,1) = urot(:)
270  vecpair(k,:,2) = vrot(:)
271  end do
272 
273  if (debug) write(logunit,'(a)')'exit '//trim(subname)
274  end subroutine getvecpair3d
275 
285  subroutine getfield2d(fname, vname, dims, field, wgts)
287  character(len=*), intent(in) :: fname, vname
288  integer, intent(in) :: dims(:)
289  real(kind=8), intent(out) :: field(:)
290  character(len=*), optional, intent(in) :: wgts
291 
292  ! local variable
293  integer :: ncid, varid
294  real(kind=8), allocatable :: a2d(:,:)
295  real(kind=8), allocatable :: atmp(:)
296  character(len=20) :: subname = 'getfield2d'
297  !----------------------------------------------------------------------------
298 
299  if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname
300 
301  allocate(a2d(dims(1),dims(2))); a2d = 0.0
302  allocate(atmp(dims(1)*dims(2))); atmp = 0.0
303 
304  call nf90_err(nf90_open(fname, nf90_nowrite, ncid), 'nf90_open: '//fname)
305  call nf90_err(nf90_inq_varid(ncid, vname, varid), &
306  'get variable ID: '//vname)
307  call nf90_err(nf90_get_var(ncid, varid, a2d), &
308  'get variable: '//vname)
309  call nf90_err(nf90_close(ncid), 'close: '//fname)
310 
311  atmp(:) = reshape(a2d, (/dims(1)*dims(2)/))
312  if(present(wgts)) then
313  call remap(trim(wgts), src_field=atmp, dst_field=field)
314  else
315  field = atmp
316  end if
317 
318  if (debug) write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname
319  end subroutine getfield2d
320 
330  subroutine getfield3d(fname, vname, dims, field, wgts)
332  character(len=*), intent(in) :: fname, vname
333  integer, intent(in) :: dims(:)
334  real(kind=8), intent(out) :: field(:,:)
335  character(len=*), optional, intent(in) :: wgts
336 
337  ! local variable
338  integer :: k, ncid, varid
339  real(kind=8), allocatable :: a3d(:,:,:)
340  real(kind=8), allocatable :: atmp(:,:)
341  character(len=20) :: subname = 'getfield3d'
342  !----------------------------------------------------------------------------
343 
344  if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname
345 
346  allocate(a3d(dims(1),dims(2),dims(3))); a3d = 0.0
347  allocate(atmp(dims(3),dims(1)*dims(2))); atmp = 0.0
348 
349  call nf90_err(nf90_open(fname, nf90_nowrite, ncid), 'nf90_open: '//fname)
350  call nf90_err(nf90_inq_varid(ncid, vname, varid), &
351  'get variable ID: '//vname)
352  call nf90_err(nf90_get_var(ncid, varid, a3d), &
353  'get variable: '//vname)
354  call nf90_err(nf90_close(ncid), 'close: '//fname)
355 
356  do k = 1,dims(3)
357  atmp(k,:) = reshape(a3d(1:dims(1),1:dims(2),k), (/dims(1)*dims(2)/))
358  end do
359  if(present(wgts)) then
360  call remap(trim(wgts), dim2=dims(3), src_field=atmp, dst_field=field)
361  else
362  field = atmp
363  end if
364 
365  if (debug) write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname
366  end subroutine getfield3d
367 
375  subroutine remap1d(fname, src_field, dst_field)
377  character(len=*), intent(in) :: fname
378  real(kind=8), intent(in) :: src_field(:)
379  real(kind=8), intent(out) :: dst_field(:)
380 
381  ! local variables
382  integer :: ncid, id
383  integer :: i,ii,jj
384  integer :: n_a, n_b, n_s
385  integer(kind=4), allocatable, dimension(:) :: col, row
386  real(kind=8), allocatable, dimension(:) :: S
387  character(len=20) :: subname = 'remap1d'
388  !----------------------------------------------------------------------------
389 
390  if (debug)write(logunit,'(a)')'enter '//trim(subname)
391 
392  ! retrieve the weights
393  call nf90_err(nf90_open(trim(fname), nf90_nowrite, ncid), 'open: '//fname)
394  call nf90_err(nf90_inq_dimid(ncid, 'n_s', id), 'get dimension Id: n_s')
395  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_s), 'get dimension: n_s' )
396  call nf90_err(nf90_inq_dimid(ncid, 'n_a', id), 'get dimension Id: n_a')
397  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_a), 'get dimension: n_a' )
398  call nf90_err(nf90_inq_dimid(ncid, 'n_b', id), 'get dimension Id: n_b')
399  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_b), 'get dimension: n_b' )
400 
401  allocate(col(1:n_s)); col = 0
402  allocate(row(1:n_s)); row = 0
403  allocate( s(1:n_s)); s = 0.0
404 
405  call nf90_err(nf90_inq_varid(ncid, 'col', id),'get variable Id: col')
406  call nf90_err(nf90_get_var(ncid, id, col),'get variable: col')
407  call nf90_err(nf90_inq_varid(ncid, 'row', id),'get variable Id: row')
408  call nf90_err(nf90_get_var(ncid, id, row),'get variable: row')
409  call nf90_err(nf90_inq_varid(ncid, 'S', id),'get variable Id: S')
410  call nf90_err(nf90_get_var(ncid, id, s),'get variable: S')
411  call nf90_err(nf90_close(ncid), 'close: '//fname)
412 
413  dst_field = 0.0
414  do i = 1,n_s
415  ii = row(i); jj = col(i)
416  dst_field(ii) = dst_field(ii) + s(i)*src_field(jj)
417  enddo
418 
419  if (debug) write(logunit,'(a)')'exit '//trim(subname)
420  end subroutine remap1d
421 
430  subroutine remap2d(fname, dim2, src_field, dst_field)
432  character(len=*), intent(in) :: fname
433  integer, intent(in) :: dim2
434  real(kind=8), intent(in) :: src_field(:,:)
435  real(kind=8), intent(out) :: dst_field(:,:)
436 
437  ! local variables
438  integer :: ncid, id
439  integer :: i,ii,jj
440  integer :: n_a, n_b, n_s
441  integer(kind=4), allocatable, dimension(:) :: col, row
442  real(kind=8), allocatable, dimension(:) :: S
443  character(len=20) :: subname = 'remap2d'
444  !----------------------------------------------------------------------------
445 
446  if (debug)write(logunit,'(a)')'enter '//trim(subname)//' weights = '//trim(fname)
447 
448  ! retrieve the weights
449  call nf90_err(nf90_open(trim(fname), nf90_nowrite, ncid), 'open: '//fname)
450  call nf90_err(nf90_inq_dimid(ncid, 'n_s', id), 'get dimension Id: n_s')
451  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_s), 'get dimension: n_s')
452  call nf90_err(nf90_inq_dimid(ncid, 'n_a', id), 'get dimension Id: n_a')
453  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_a), 'get dimension: n_a')
454  call nf90_err(nf90_inq_dimid(ncid, 'n_b', id), 'get dimension Id: n_b')
455  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_b), 'get dimension: n_b')
456 
457  allocate(col(1:n_s)); col = 0
458  allocate(row(1:n_s)); row = 0
459  allocate( s(1:n_s)); s = 0.0
460 
461  call nf90_err(nf90_inq_varid(ncid, 'col', id),'get variable Id: col')
462  call nf90_err(nf90_get_var(ncid, id, col),'get variable: col')
463  call nf90_err(nf90_inq_varid(ncid, 'row', id),'get variable Id: row')
464  call nf90_err(nf90_get_var(ncid, id, row),'get variable: row')
465  call nf90_err(nf90_inq_varid(ncid, 'S', id),'get variable Id: S')
466  call nf90_err(nf90_get_var(ncid, id, s),'get variable: S')
467  call nf90_err(nf90_close(ncid), 'close: '//fname)
468 
469  dst_field = 0.0
470  do i = 1,n_s
471  ii = row(i); jj = col(i)
472  dst_field(:,ii) = dst_field(:,ii) + s(i)*src_field(:,jj)
473  enddo
474 
475  if (debug) write(logunit,'(a)')'exit '//trim(subname)
476  end subroutine remap2d
477 
487  subroutine remap3d(fname, nk, nflds, src_field, dst_field)
489  character(len=*), intent(in) :: fname
490  integer, intent(in) :: nk, nflds
491  real(kind=8), intent(in) :: src_field(:,:,:)
492  real(kind=8), intent(out) :: dst_field(:,:,:)
493 
494  ! local variables
495  integer :: ncid, id
496  integer :: i,ii,jj
497  integer :: n_a, n_b, n_s
498  integer(kind=4), allocatable, dimension(:) :: col, row
499  real(kind=8), allocatable, dimension(:) :: S
500  character(len=20) :: subname = 'remap3d'
501  !----------------------------------------------------------------------------
502 
503  if (debug)write(logunit,'(a)')'enter '//trim(subname)//' weights = '//trim(fname)
504 
505  ! retrieve the weights
506  call nf90_err(nf90_open(trim(fname), nf90_nowrite, ncid), 'open: '//fname)
507  call nf90_err(nf90_inq_dimid(ncid, 'n_s', id), 'get dimension Id: n_s')
508  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_s), 'get dimension: n_s')
509  call nf90_err(nf90_inq_dimid(ncid, 'n_a', id), 'get dimension Id: n_a')
510  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_a), 'get dimension: n_a')
511  call nf90_err(nf90_inq_dimid(ncid, 'n_b', id), 'get dimension Id: n_b')
512  call nf90_err(nf90_inquire_dimension(ncid, id, len=n_b), 'get dimension: n_b')
513 
514  allocate(col(1:n_s)); col = 0
515  allocate(row(1:n_s)); row = 0
516  allocate( s(1:n_s)); s = 0.0
517 
518  call nf90_err(nf90_inq_varid(ncid, 'col', id),'get variable Id: col')
519  call nf90_err(nf90_get_var(ncid, id, col),'get variable: col')
520  call nf90_err(nf90_inq_varid(ncid, 'row', id),'get variable Id: row')
521  call nf90_err(nf90_get_var(ncid, id, row),'get variable: row')
522  call nf90_err(nf90_inq_varid(ncid, 'S', id),'get variable Id: S')
523  call nf90_err(nf90_get_var(ncid, id, s),'get variable: S')
524  call nf90_err(nf90_close(ncid), 'close: '//fname)
525 
526  dst_field = 0.0
527  do i = 1,n_s
528  ii = row(i); jj = col(i)
529  dst_field(:,:,ii) = dst_field(:,:,ii) + s(i)*src_field(:,:,jj)
530  enddo
531 
532  if (debug) write(logunit,'(a)')'exit '//trim(subname)
533  end subroutine remap3d
534 
544  subroutine dumpnc2d(fname, vname, dims, nflds, field)
546  character(len=*), intent(in) :: fname, vname
547  integer, intent(in) :: dims(:)
548  integer, intent(in) :: nflds
549  real(kind=8), intent(in) :: field(:,:)
550 
551  ! local variable
552  integer :: n, ncid, varid, idimid, jdimid, fdimid
553  real(kind=8), allocatable :: a3d(:,:,:)
554  character(len=20) :: subname = 'dumpnc2d'
555  !----------------------------------------------------------------------------
556 
557  if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname
558  allocate(a3d(dims(1),dims(2),nflds)); a3d = 0.0
559 
560  call nf90_err(nf90_create(trim(fname), nf90_clobber, ncid), 'create: '//fname)
561  call nf90_err(nf90_def_dim(ncid, 'nx', dims(1), idimid), 'define dimension: nx')
562  call nf90_err(nf90_def_dim(ncid, 'ny', dims(2), jdimid), 'define dimension: ny')
563  call nf90_err(nf90_def_dim(ncid, 'nf', nflds, fdimid), 'define dimension: nf')
564  call nf90_err(nf90_def_var(ncid, vname, nf90_double, (/idimid,jdimid,fdimid/), varid), &
565  'define variable: '//vname)
566  call nf90_err(nf90_enddef(ncid), 'nf90_enddef: '//fname)
567 
568  do n = 1,nflds
569  a3d(:,:,n) = reshape(field(n,1:dims(1)*dims(2)), (/dims(1),dims(2)/))
570  end do
571  call nf90_err(nf90_put_var(ncid, varid, a3d), 'put variable: '//vname)
572  call nf90_err(nf90_close(ncid), 'close: '//fname)
573 
574  if (debug)write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname
575  end subroutine dumpnc2d
576 
587  subroutine dumpnc3d(fname, vname, dims, nk, nflds, field)
589  character(len=*), intent(in) :: fname, vname
590  integer, intent(in) :: dims(:)
591  integer, intent(in) :: nk, nflds
592  real(kind=8), intent(in) :: field(:,:,:)
593 
594  ! local variable
595  integer :: n, k, ncid, varid, idimid, jdimid, kdimid, fdimid
596  real(kind=8), allocatable :: a4d(:,:,:,:)
597  character(len=20) :: subname = 'dumpnc3d'
598  !----------------------------------------------------------------------------
599 
600  if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname
601  allocate(a4d(dims(1),dims(2),dims(3),nflds)); a4d = 0.0
602 
603  call nf90_err(nf90_create(trim(fname), nf90_clobber, ncid), 'create: '//fname)
604  call nf90_err(nf90_def_dim(ncid, 'nx', dims(1), idimid), 'define dimension: nx')
605  call nf90_err(nf90_def_dim(ncid, 'ny', dims(2), jdimid), 'define dimension: ny')
606  call nf90_err(nf90_def_dim(ncid, 'nk', dims(3), kdimid), 'define dimension: nk')
607  call nf90_err(nf90_def_dim(ncid, 'nf', nflds, fdimid), 'define dimension: nf')
608  call nf90_err(nf90_def_var(ncid, vname, nf90_double, (/idimid,jdimid,kdimid,fdimid/), varid), &
609  'define variable: '//vname)
610  call nf90_err(nf90_enddef(ncid), 'nf90_enddef: '//fname)
611 
612  do n = 1,nflds
613  do k = 1,nk
614  a4d(:,:,k,n) = reshape(field(n,k,1:dims(1)*dims(2)), (/dims(1),dims(2)/))
615  end do
616  end do
617  call nf90_err(nf90_put_var(ncid, varid, a4d), 'put variable: '//vname)
618  call nf90_err(nf90_close(ncid), 'close: '//fname)
619 
620  if (debug)write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname
621  end subroutine dumpnc3d
622 
631  subroutine dumpnc3dk(fname, vname, dims, field)
633  character(len=*), intent(in) :: fname, vname
634  integer, intent(in) :: dims(:)
635  real(kind=8), intent(in) :: field(:,:)
636 
637  ! local variable
638  integer :: k, ncid, varid, idimid, jdimid, kdimid
639  real(kind=8), allocatable :: a3d(:,:,:)
640  character(len=20) :: subname = 'dumpnc3dk'
641  !----------------------------------------------------------------------------
642 
643  if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname
644  allocate(a3d(dims(1),dims(2),dims(3))); a3d = 0.0
645 
646  call nf90_err(nf90_create(trim(fname), nf90_clobber, ncid), 'nf90_create: '//fname)
647  call nf90_err(nf90_def_dim(ncid, 'nx', dims(1), idimid), 'define dimension: nx')
648  call nf90_err(nf90_def_dim(ncid, 'ny', dims(2), jdimid), 'define dimension: ny')
649  call nf90_err(nf90_def_dim(ncid, 'nk', dims(3), kdimid), 'define dimension: nk')
650  call nf90_err(nf90_def_var(ncid, vname, nf90_double, (/idimid,jdimid,kdimid/), varid), &
651  'define variable: '//vname)
652  call nf90_err(nf90_enddef(ncid), 'nf90_enddef: '//fname)
653 
654  do k = 1,dims(3)
655  a3d(:,:,k) = reshape(field(k,1:dims(1)*dims(2)), (/dims(1),dims(2)/))
656  end do
657  call nf90_err(nf90_put_var(ncid, varid, a3d), 'put variable: '//vname)
658  call nf90_err(nf90_close(ncid), 'close: '//fname)
659 
660  if (debug)write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname
661 
662  end subroutine dumpnc3dk
663 
672  subroutine dumpnc1d(fname, vname, dims, field)
674  character(len=*), intent(in) :: fname, vname
675  integer, intent(in) :: dims(:)
676  real(kind=8), intent(in) :: field(:)
677 
678  ! local variable
679  integer :: ncid, varid, idimid, jdimid
680  real(kind=8), allocatable :: a2d(:,:)
681  character(len=20) :: subname = 'dumpnc1d'
682  !----------------------------------------------------------------------------
683 
684  if (debug)write(logunit,'(a)')'enter '//trim(subname)//' variable '//vname
685  allocate(a2d(dims(1),dims(2))); a2d = 0.0
686 
687  call nf90_err(nf90_create(trim(fname), nf90_clobber, ncid), 'nf90_create: '//fname)
688  call nf90_err(nf90_def_dim(ncid, 'nx', dims(1), idimid), 'define dimension: nx')
689  call nf90_err(nf90_def_dim(ncid, 'ny', dims(2), jdimid), 'define dimension: ny')
690  call nf90_err(nf90_def_var(ncid, vname, nf90_double, (/idimid,jdimid/), varid), &
691  'define variable: '//vname)
692  call nf90_err(nf90_enddef(ncid), 'nf90_enddef: '//fname)
693 
694  a2d(:,:) = reshape(field(1:dims(1)*dims(2)), (/dims(1),dims(2)/))
695  call nf90_err(nf90_put_var(ncid, varid, a2d), 'put variable: '//vname)
696  call nf90_err(nf90_close(ncid), 'close: '//fname)
697 
698  if (debug)write(logunit,'(a)')'exit '//trim(subname)//' variable '//vname
699 
700  end subroutine dumpnc1d
701 
708  subroutine nf90_err(ierr, string)
710  integer , intent(in) :: ierr
711  character(len=*), intent(in) :: string
712  !----------------------------------------------------------------------------
713 
714  if (ierr /= nf90_noerr) then
715  write(0, '(a)') 'FATAL ERROR: ' // trim(string)// ' : ' // trim(nf90_strerror(ierr))
716  ! This fails on WCOSS2 with Intel 19 compiler. See
717  ! https://community.intel.com/t5/Intel-Fortran-Compiler/STOP-and-ERROR-STOP-with-variable-stop-codes/m-p/1182521#M149254
718  ! When WCOSS2 moves to Intel 2020+, uncomment the next line and remove stop 99
719  !stop ierr
720  stop 99
721  end if
722  end subroutine nf90_err
723 end module utils_mod