ocnice_prep  1.13.0
•All Data Structures Files Functions Variables Pages
utils_esmf_mod.F90
Go to the documentation of this file.
1 
7 module utils_esmf_mod
8 
9  use esmf
10  use netcdf
11  use init_mod , only : do_ocnprep, debug, logunit, vardefs, fsrc, fdst, ftype, nxr, nyr
12  use utils_mod , only : dumpnc, remap
13  use arrays_mod, only : maskspval, mask3d
14 
15  implicit none
16 
17  private
18 
19  type(esmf_routehandle) :: rh
20  type(esmf_dynamicmask) :: dynamiclevmask
21  type(esmf_mesh) :: meshsrc
22  type(esmf_mesh) :: meshdst
23 
24  integer :: srctermprocessing = 0
25 
26  interface remaprh
27  module procedure remaprh1d
28  module procedure remaprh2d
29  module procedure remaprh1ddyn
30  module procedure remaprh2ddyn
31  end interface remaprh
32 
33  interface rotremap
34  module procedure rotremap2d
35  module procedure rotremap3d
36  end interface rotremap
37 
38  public :: createrh
39  public :: remaprh
40  public :: rotremap
41  public :: chkerr
42 
43  character(len=*), parameter :: u_file_u = __file__
44 contains
52  subroutine createrh(srcmeshfile,dstmeshfile,rc)
53 
54  character(len=*), intent(in) :: srcmeshfile
55  character(len=*), intent(in) :: dstmeshfile
56  integer, intent(out) :: rc
57 
58  ! local variables
59  type(esmf_field) :: fldsrc
60  type(esmf_field) :: flddst
61  type(esmf_regridmethod_flag) :: regridmethod
62  type(esmf_extrapmethod_flag) :: extrapmethod
63  type(esmf_field) :: dststatusfield
64  integer, pointer :: dststatus(:)
65  character(len=20) :: subname = 'remapRH1d'
66  !----------------------------------------------------------------------------
67 
68  if (debug)write(logunit,'(a)')'enter '//trim(subname)
69  rc = esmf_success
70 
71  ! use nstod for ice to maintain thermodynamically consistent values
72  ! w/in layers and categories
73  if (do_ocnprep) then
74  regridmethod = esmf_regridmethod_bilinear
75  extrapmethod = esmf_extrapmethod_nearest_idavg
76  else
77  regridmethod = esmf_regridmethod_nearest_stod
78  extrapmethod = esmf_extrapmethod_nearest_stod
79  end if
80 
81  meshsrc = esmf_meshcreate(filename=trim(srcmeshfile), fileformat=esmf_fileformat_esmfmesh, rc=rc)
82  if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
83  fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, name='mshsrc', meshloc=esmf_meshloc_element, rc=rc)
84  if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
85 
86  meshdst = esmf_meshcreate(filename=trim(dstmeshfile), fileformat=esmf_fileformat_esmfmesh, rc=rc)
87  if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
88  flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, name='mshdst', meshloc=esmf_meshloc_element, rc=rc)
89  if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
90 
91  dststatusfield = esmf_fieldcreate(meshdst, esmf_typekind_i4, meshloc=esmf_meshloc_element, rc=rc)
92  if (chkerr(rc,__line__,u_file_u)) return
93  call esmf_fieldget(dststatusfield, farrayptr=dststatus, rc=rc)
94  if (chkerr(rc,__line__,u_file_u)) return
95 
96  call esmf_fieldregridstore(fldsrc, flddst, routehandle=rh, &
97  srcmaskvalues=(/0/), dstmaskvalues=(/0/), &
98  regridmethod=regridmethod, &
99  extrapmethod=extrapmethod, &
100  polemethod=esmf_polemethod_allavg, &
101  ignoredegenerate=.true., &
102  srctermprocessing=srctermprocessing, &
103  dststatusfield=dststatusfield, &
104  unmappedaction=esmf_unmappedaction_ignore, rc=rc)
105  if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
106 
107  ! Create a dynamic mask object
108  call esmf_dynamicmasksetr8r8r8(dynamiclevmask, dynamicsrcmaskvalue=maskspval, &
109  dynamicmaskroutine=dynlevmaskproc, rc=rc)
110  if (chkerr(rc,__line__,u_file_u)) return
111 
112  if (debug) then
113  call dumpnc(trim(ftype)//'.'//trim(fdst)//'.dststat.nc', 'dststat', &
114  dims=(/nxr,nyr/), field=real(dststatus,8))
115  end if
116 
117  if (debug)write(logunit,'(a)')'exit '//trim(subname)
118  end subroutine createrh
119 
128  subroutine remaprh1d(kk,src_field,dst_field,rc)
129 
130  integer, intent(in) :: kk
131  real(kind=8), intent(in) :: src_field(:)
132  real(kind=8), intent(out) :: dst_field(:)
133  integer, intent(out) :: rc
134 
135  ! local variables
136  type(ESMF_Field) :: fldsrc
137  type(ESMF_Field) :: flddst
138  real(kind=8), pointer :: srcptr(:), dstptr(:)
139  character(len=20) :: subname = 'remapRH1d'
140  !----------------------------------------------------------------------------
141 
142  if (debug)write(logunit,'(a,i5)')'enter '//trim(subname)//' ',kk
143  rc = esmf_success
144 
145  fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, meshloc=esmf_meshloc_element,rc=rc)
146  if (chkerr(rc,__line__,u_file_u)) return
147  flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, meshloc=esmf_meshloc_element,rc=rc)
148  if (chkerr(rc,__line__,u_file_u)) return
149 
150  call esmf_fieldget(fldsrc, farrayptr=srcptr, rc=rc)
151  if (chkerr(rc,__line__,u_file_u)) return
152  call esmf_fieldget(flddst, farrayptr=dstptr, rc=rc)
153  if (chkerr(rc,__line__,u_file_u)) return
154 
155  srcptr = src_field
156  if (esmf_routehandleiscreated(rh,rc=rc)) then
157  call esmf_fieldregrid(fldsrc, flddst, routehandle=rh, rc=rc)
158  if (chkerr(rc,__line__,u_file_u)) return
159  else
160  call esmf_logwrite(trim(subname)//": RH not created ", esmf_logmsg_info)
161  rc=esmf_failure
162  end if
163  dst_field = dstptr
164 
165  if (debug)write(logunit,'(a,i5)')'exit '//trim(subname)//' ',kk
166  end subroutine remaprh1d
167 
175  subroutine remaprh2d(src_field,dst_field,rc)
176 
177  real(kind=8), intent(in) :: src_field(:,:)
178  real(kind=8), intent(out) :: dst_field(:,:)
179  integer, intent(out) :: rc
180 
181  ! local variables
182  type(ESMF_Field) :: fldsrc
183  type(ESMF_Field) :: flddst
184  real(kind=8), pointer :: srcptr(:,:), dstptr(:,:)
185  character(len=20) :: subname = 'remapRH2d'
186  !----------------------------------------------------------------------------
187 
188  if (debug)write(logunit,'(a)')'enter '//trim(subname)
189  rc = esmf_success
190 
191  fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, meshloc=esmf_meshloc_element, &
192  ungriddedlbound=(/1/), ungriddedubound=(/size(src_field,1)/), &
193  gridtofieldmap=(/2/), rc=rc)
194  if (chkerr(rc,__line__,u_file_u)) return
195  flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, meshloc=esmf_meshloc_element, &
196  ungriddedlbound=(/1/), ungriddedubound=(/size(dst_field,1)/), &
197  gridtofieldmap=(/2/), rc=rc)
198  if (chkerr(rc,__line__,u_file_u)) return
199 
200  call esmf_fieldget(fldsrc, farrayptr=srcptr, rc=rc)
201  if (chkerr(rc,__line__,u_file_u)) return
202  call esmf_fieldget(flddst, farrayptr=dstptr, rc=rc)
203  if (chkerr(rc,__line__,u_file_u)) return
204 
205  srcptr = src_field
206  if (esmf_routehandleiscreated(rh,rc=rc)) then
207  call esmf_fieldregrid(fldsrc, flddst, routehandle=rh, rc=rc)
208  if (chkerr(rc,__line__,u_file_u)) return
209  else
210  call esmf_logwrite(trim(subname)//": RH not created ", esmf_logmsg_info)
211  rc=esmf_failure
212  end if
213  dst_field = dstptr
214 
215  if (debug)write(logunit,'(a)')'exit '//trim(subname)
216  end subroutine remaprh2d
217 
227  subroutine remaprh1ddyn(kk,src_field,dst_field,hmask,rc)
228 
229  !nflds,nlen
230  integer, intent(in) :: kk
231  real(kind=8), intent(in) :: src_field(:)
232  real(kind=8), intent(in) :: hmask(:)
233  real(kind=8), intent(out) :: dst_field(:)
234  integer, intent(out) :: rc
235 
236  ! local variables
237  type(ESMF_Field) :: fldsrc
238  type(ESMF_Field) :: flddst
239  real(kind=8), pointer :: srcptr(:), dstptr(:)
240  character(len=20) :: subname = 'remapRH1ddyn'
241  !----------------------------------------------------------------------------
242 
243  if (debug)write(logunit,'(a,i5)')'enter '//trim(subname)//' ',kk
244  rc = esmf_success
245 
246  fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, meshloc=esmf_meshloc_element)
247  if (chkerr(rc,__line__,u_file_u)) return
248  flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, meshloc=esmf_meshloc_element)
249  if (chkerr(rc,__line__,u_file_u)) return
250 
251  call esmf_fieldget(fldsrc, farrayptr=srcptr, rc=rc)
252  if (chkerr(rc,__line__,u_file_u)) return
253  call esmf_fieldget(flddst, farrayptr=dstptr, rc=rc)
254  if (chkerr(rc,__line__,u_file_u)) return
255 
256  srcptr = src_field
257  where(hmask .eq. maskspval)srcptr = maskspval
258 
259  if (esmf_routehandleiscreated(rh,rc=rc)) then
260  call esmf_fieldregrid(fldsrc, flddst, routehandle=rh, dynamicmask=dynamiclevmask, rc=rc)
261  if (chkerr(rc,__line__,u_file_u)) return
262  else
263  call esmf_logwrite(trim(subname)//": RH not created ", esmf_logmsg_info)
264  rc=esmf_failure
265  end if
266  dst_field = dstptr
267 
268  if (debug)write(logunit,'(a,i5)')'exit '//trim(subname)//' ',kk
269  end subroutine remaprh1ddyn
270 
280  subroutine remaprh2ddyn(kk,src_field,dst_field,hmask,rc)
281 
282  integer, intent(in) :: kk
283  real(kind=8), intent(in) :: src_field(:,:)
284  real(kind=8), intent(in) :: hmask(:)
285  real(kind=8), intent(out) :: dst_field(:,:)
286  integer, intent(out) :: rc
287 
288  ! local variables
289  type(ESMF_Field) :: fldsrc
290  type(ESMF_Field) :: flddst
291  integer :: i,n
292  real(kind=8), pointer :: srcptr(:,:), dstptr(:,:)
293  character(len=20) :: subname = 'remapRH2ddyn'
294  !----------------------------------------------------------------------------
295 
296  if (debug)write(logunit,'(a,i5)')'enter '//trim(subname)//' ',kk
297  rc = esmf_success
298 
299  fldsrc = esmf_fieldcreate(meshsrc, esmf_typekind_r8, meshloc=esmf_meshloc_element, &
300  ungriddedlbound=(/1/), ungriddedubound=(/size(src_field,1)/), &
301  gridtofieldmap=(/2/), rc=rc)
302  if (chkerr(rc,__line__,u_file_u)) return
303  flddst = esmf_fieldcreate(meshdst, esmf_typekind_r8, meshloc=esmf_meshloc_element, &
304  ungriddedlbound=(/1/), ungriddedubound=(/size(dst_field,1)/), &
305  gridtofieldmap=(/2/), rc=rc)
306  if (chkerr(rc,__line__,u_file_u)) return
307 
308  call esmf_fieldget(fldsrc, farrayptr=srcptr, rc=rc)
309  if (chkerr(rc,__line__,u_file_u)) return
310  call esmf_fieldget(flddst, farrayptr=dstptr, rc=rc)
311  if (chkerr(rc,__line__,u_file_u)) return
312 
313  srcptr = src_field
314  do n = 1,ubound(src_field,2)
315  do i = 1,ubound(src_field,1)
316  if(hmask(n) .eq. maskspval)srcptr(i,n) = maskspval
317  end do
318  end do
319 
320  if (esmf_routehandleiscreated(rh,rc=rc)) then
321  call esmf_fieldregrid(fldsrc, flddst, routehandle=rh, dynamicmask=dynamiclevmask, rc=rc)
322  if (chkerr(rc,__line__,u_file_u)) return
323  else
324  call esmf_logwrite(trim(subname)//": RH not created ", esmf_logmsg_info)
325  rc=esmf_failure
326  end if
327  dst_field = dstptr
328 
329  if (debug)write(logunit,'(a,i5)')'exit '//trim(subname)//' ',kk
330  end subroutine remaprh2ddyn
331 
343  subroutine rotremap2d(wdir, vars, cosrot, sinrot, dims, nflds, fields)
344 
345  character(len=*), intent(in) :: wdir
346  real(kind=8), intent(in) :: cosrot(:),sinrot(:)
347  type(vardefs), intent(in) :: vars(:)
348  integer, intent(in) :: dims(:)
349  integer, intent(in) :: nflds
350  real(kind=8), intent(inout) :: fields(:,:)
351 
352  integer :: n, idx1, idx2
353  real(kind=8), allocatable, dimension(:) :: urot, vrot
354  character(len=10) :: vgrid1, vgrid2
355  character(len=240) :: wgtsfile
356  character(len=20) :: subname = 'rotremap2d'
357  !----------------------------------------------------------------------------
358 
359  if (debug)write(logunit,'(a)')'enter '//trim(subname)
360 
361  idx1 = 0; idx2 = 0
362  do n = 1,nflds
363  if (len_trim(vars(n)%var_pair) > 0 .and. idx1 .eq. 0) then
364  idx1 = n
365  idx2 = n+1
366  end if
367  end do
368  if (idx1 .eq. 0)return
369 
370  vgrid1 = vars(idx1)%var_grid(1:2)
371  vgrid2 = vars(idx1)%var_pair_grid(1:2)
372 
373  allocate(urot(1:dims(1)*dims(2))); urot = 0.0
374  allocate(vrot(1:dims(1)*dims(2))); vrot = 0.0
375  urot(:) = fields(idx1,:)*cosrot(:) - fields(idx2,:)*sinrot(:)
376  vrot(:) = fields(idx2,:)*cosrot(:) + fields(idx1,:)*sinrot(:)
377 
378  wgtsfile = trim(wdir)//'tripole.'//trim(fdst)//'.Ct.to.'//trim(vgrid1)//'.bilinear.nc'
379  call remap(trim(wgtsfile), urot, fields(idx1,:))
380  if (debug) write(logunit,'(a)')'use '//trim(wgtsfile)//' to restagger from Ct to '//trim(vgrid1)
381 
382  wgtsfile = trim(wdir)//'tripole.'//trim(fdst)//'.Ct.to.'//trim(vgrid2)//'.bilinear.nc'
383  call remap(trim(wgtsfile), vrot, fields(idx2,:))
384  if (debug) write(logunit,'(a)')'use '//trim(wgtsfile)//' to restagger from Ct to '//trim(vgrid2)
385 
386  if (debug)write(logunit,'(a)')'exit '//trim(subname)
387  end subroutine rotremap2d
388 
400  subroutine rotremap3d(wdir, vars, cosrot, sinrot, dims, nflds, fields)
401 
402  character(len=*), intent(in) :: wdir
403  real(kind=8), intent(in) :: cosrot(:),sinrot(:)
404  type(vardefs), intent(in) :: vars(:)
405  integer, intent(in) :: dims(:)
406  integer, intent(in) :: nflds
407  real(kind=8), intent(inout) :: fields(:,:,:)
408 
409  integer :: k, n, idx1, idx2
410  real(kind=8), allocatable, dimension(:) :: urot, vrot
411  character(len=10) :: vgrid1, vgrid2
412  character(len=240) :: wgtsfile
413  character(len=20) :: subname = 'rotremap3d'
414  !----------------------------------------------------------------------------
415 
416  if (debug)write(logunit,'(a)')'enter '//trim(subname)
417 
418  idx1 = 0; idx2 = 0
419  do n = 1,nflds
420  if (len_trim(vars(n)%var_pair) > 0 .and. idx1 .eq. 0) then
421  idx1 = n
422  idx2 = n+1
423  end if
424  end do
425  if (idx1 .eq. 0)return
426 
427  vgrid1 = vars(idx1)%var_grid(1:2)
428  vgrid2 = vars(idx1)%var_pair_grid(1:2)
429  if (debug) write(logunit,'(a)')'restagger from Ct to '//trim(vgrid1)//' and '//trim(vgrid2)
430 
431  allocate(urot(1:dims(1)*dims(2))); urot = 0.0
432  allocate(vrot(1:dims(1)*dims(2))); vrot = 0.0
433  do k = 1,dims(3)
434  urot(:) = fields(idx1,k,:)*cosrot(:) - fields(idx2,k,:)*sinrot(:)
435  vrot(:) = fields(idx2,k,:)*cosrot(:) + fields(idx1,k,:)*sinrot(:)
436  wgtsfile = trim(wdir)//'tripole.'//trim(fdst)//'.Ct.to.'//trim(vgrid1)//'.bilinear.nc'
437  call remap(trim(wgtsfile), urot, fields(idx1,k,:))
438  wgtsfile = trim(wdir)//'tripole.'//trim(fdst)//'.Ct.to.'//trim(vgrid2)//'.bilinear.nc'
439  call remap(trim(wgtsfile), vrot, fields(idx2,k,:))
440  end do
441 
442  if (debug)write(logunit,'(a)')'exit '//trim(subname)
443  end subroutine rotremap3d
444 
453  subroutine dynlevmaskproc(dynamicMaskList, dynamicSrcMaskValue, dynamicDstMaskValue, rc)
454 
455  ! input/output arguments
456  type(ESMF_DynamicMaskElementR8R8R8) , pointer :: dynamicMaskList(:)
457  real(ESMF_KIND_R8), intent(in), optional :: dynamicSrcMaskValue
458  real(ESMF_KIND_R8), intent(in), optional :: dynamicDstMaskValue
459  integer , intent(out) :: rc
460 
461  ! local variables
462  integer :: i, j
463  real(ESMF_KIND_R8) :: renorm
464  character(len=20) :: subname = 'dynLevMaskProc'
465  !----------------------------------------------------------------------------
466 
467  if (debug)write(logunit,'(a)')'enter '//trim(subname)
468  rc = esmf_success
469 
470  if (associated(dynamicmasklist)) then
471  do i=1, size(dynamicmasklist)
472  dynamicmasklist(i)%dstElement = 0.0 ! set to zero
473  renorm = 0.0 ! reset
474  do j = 1, size(dynamicmasklist(i)%factor)
475  if (dynamicsrcmaskvalue /= dynamicmasklist(i)%srcElement(j)) then
476  dynamicmasklist(i)%dstElement = dynamicmasklist(i)%dstElement + &
477  (dynamicmasklist(i)%factor(j) * dynamicmasklist(i)%srcElement(j))
478  renorm = renorm + dynamicmasklist(i)%factor(j)
479  endif
480  enddo
481  if (renorm > 0.0) then
482  dynamicmasklist(i)%dstElement = dynamicmasklist(i)%dstElement / renorm
483  else if (present(dynamicsrcmaskvalue)) then
484  dynamicmasklist(i)%dstElement = dynamicsrcmaskvalue
485  else
486  rc = esmf_rc_arg_bad ! error detected
487  return
488  endif
489  enddo
490  endif
491  if (debug)write(logunit,'(a)')'exit '//trim(subname)
492 
493  end subroutine dynlevmaskproc
494 
503  logical function chkerr(rc, line, file)
504  integer, intent(in) :: rc
505  integer, intent(in) :: line
506  character(len=*), intent(in) :: file
507  integer :: lrc
508  !----------------------------------------------------------------------------
509  chkerr = .false.
510  lrc = rc
511  if (esmf_logfounderror(rctocheck=lrc, msg=esmf_logerr_passthru, line=line, file=file)) then
512  chkerr = .true.
513  endif
514  end function chkerr
515 end module utils_esmf_mod
utils_esmf_mod::rotremap
Definition: utils_esmf_mod.F90:33
utils_esmf_mod::remaprh
Definition: utils_esmf_mod.F90:26
utils_mod::dumpnc
Definition: utils_mod.F90:38
utils_mod::remap
Definition: utils_mod.F90:27
init_mod::vardefs
Definition: init_mod.F90:16