ocnice_prep 1.14.0
Loading...
Searching...
No Matches
ocniceprep.F90
Go to the documentation of this file.
1
6
18program ocniceprep
19
20 use esmf
21 use netcdf
22 use init_mod , only : nxt, nyt, nlevs, nxr, nyr, outvars, readnml, readcsv
23 use init_mod , only : wgtsdir, griddir, ftype, fsrc, fdst, input_file, maskvar
24 use init_mod , only : do_ocnprep, debug, logunit
25 use arrays_mod , only : b2d, c2d, b3d, rgb2d, rgb3d, rgc2d, setup_packing
26 use arrays_mod , only : nbilin2d, nbilin3d, nconsd2d, bilin2d, bilin3d, consd2d
27 use arrays_mod , only : mask3d, hmin, maskspval, eta
28 use utils_mod , only : getfield, packarrays, remap, dumpnc, nf90_err
29 use utils_mod , only : zero_out_land_ice, zero_out_phantom_ice
30 use utils_esmf_mod , only : createrh, remaprh, chkerr, rotremap
31 use restarts_mod , only : setup_icerestart, setup_ocnrestart
32 use ocncalc_mod , only : calc_eta, vfill
33
34 implicit none
35
36 type(esmf_vm) :: vm
37 character(len=160) :: gridfile
38 character(len=160) :: wgtsfile
39 character(len=160) :: fout
40
41 real(kind=8), allocatable, dimension(:) :: angsrc
42 real(kind=8), allocatable, dimension(:) :: angdst
43
44 real(kind=8), allocatable, dimension(:) :: bathysrc
45 real(kind=8), allocatable, dimension(:) :: bathydst
46
47 ! work arrays for output netcdf
48 real(kind=8), allocatable, dimension(:,:) :: out2d
49 real(kind=8), allocatable, dimension(:,:,:) :: out3d
50
51 ! ice mask array for QC of ice files
52 integer(kind=4), allocatable, dimension(:,:) :: a2d
53 integer(kind=4), allocatable, dimension(:) :: kmt
54 integer :: idx1, idx2, idx3, idx4
55
56 character(len=120) :: errmsg
57 character(len=120) :: meshfsrc, meshfdst
58 integer :: nvalid, icnt
59 integer :: k,n,nn,rc,ncid,varid
60 integer :: localpet, npet
61 character(len=20) :: vname
62
63 character(len=*), parameter :: u_file_u = __file__
64
65 ! -----------------------------------------------------------------------------
66 ! initialize ESMF
67 ! -----------------------------------------------------------------------------
68
69 call esmf_vmgetglobal(vm, rc=rc)
70 call esmf_initialize(vm=vm, logkindflag=esmf_logkind_multi, rc=rc)
71 call esmf_vmget(vm, localpet=localpet, pecount=npet, rc=rc)
72 if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
73 line=__line__, file=__file__)) call esmf_finalize(endflag=esmf_end_abort)
74
75 ! -----------------------------------------------------------------------------
76 ! read the nml file and a file containing the list of variables to be remapped
77 ! -----------------------------------------------------------------------------
78
79 call readnml('ocniceprep.nml',errmsg,rc)
80 if (rc /= 0) then
81 write(0,'(a)')trim(errmsg)
82 stop
83 else
84 write(logunit,'(a)')trim(errmsg)
85 end if
86
87 call readcsv(trim(ftype)//'.csv',errmsg,rc,nvalid)
88 if (rc /= 0) then
89 write(0,'(a)')trim(errmsg)
90 stop
91 else
92 write(logunit,'(a)')trim(errmsg)
93 end if
94
95 if (npet /= 1) then
96 write(logunit, '(a)')'More than one task specified; Aborting '
97 call esmf_finalize(endflag=esmf_end_abort)
98 end if
99
100 ! -----------------------------------------------------------------------------
101 ! create a regrid RH from source to destination
102 ! -----------------------------------------------------------------------------
103
104 meshfsrc = trim(griddir)//fsrc(3:5)//'/'//'mesh.'//trim(fsrc)//'.nc'
105 meshfdst = trim(griddir)//fdst(3:5)//'/'//'mesh.'//trim(fdst)//'.nc'
106 write(logunit,'(a)')'mesh src: '//trim(meshfsrc)
107 write(logunit,'(a)')'mesh dst: '//trim(meshfdst)
108 call createrh(trim(meshfsrc),trim(meshfdst),rc=rc)
109 if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
110
111 ! -----------------------------------------------------------------------------
112 ! read the master grid file and obtain the rotation angle on the source and
113 ! destination grids
114 ! -----------------------------------------------------------------------------
115
116 allocate(angsrc(nxt*nyt)); angsrc = 0.0
117 allocate(angdst(nxr*nyr)); angdst = 0.0
118 allocate(bathysrc(nxt*nyt)); bathysrc = 0.0
119 allocate(bathydst(nxr*nyr)); bathydst = 0.0
120
121 gridfile = trim(griddir)//fsrc(3:5)//'/'//'tripole.'//trim(fsrc)//'.nc'
122 call nf90_err(nf90_open(trim(gridfile), nf90_nowrite, ncid), &
123 'open: '//trim(gridfile))
124 call getfield(trim(gridfile), 'anglet', dims=(/nxt,nyt/), field=angsrc)
125 call getfield(trim(gridfile), 'depth', dims=(/nxt,nyt/), field=bathysrc)
126 call nf90_err(nf90_close(ncid), 'close: '//trim(gridfile))
127
128 gridfile = trim(griddir)//fdst(3:5)//'/'//'tripole.'//trim(fdst)//'.nc'
129 call nf90_err(nf90_open(trim(gridfile), nf90_nowrite, ncid), &
130 'open: '//trim(gridfile))
131 call getfield(trim(gridfile), 'anglet', dims=(/nxr,nyr/), field=angdst)
132 call getfield(trim(gridfile), 'depth', dims=(/nxr,nyr/), field=bathydst)
133 call nf90_err(nf90_close(ncid), 'close: '//trim(gridfile))
134
135 if (.not. do_ocnprep) then
136 ! -----------------------------------------------------------------------------
137 ! obtain the land mask for the ice model on the source grid to QC the ice fields
138 ! -----------------------------------------------------------------------------
139
140 allocate(a2d(nxt,nyt)); a2d = 0
141 allocate(kmt(nxt*nyt)); kmt = 0
142
143 ! obtain the land mask for the source grid
144 vname = 'wet'
145 gridfile = trim(griddir)//fsrc(3:5)//'/'//'tripole.'//trim(fsrc)//'.nc'
146 call nf90_err(nf90_open(trim(gridfile), nf90_nowrite, ncid), &
147 'open: '//trim(gridfile))
148 call nf90_err(nf90_inq_varid(ncid, trim(vname), varid), 'get variable ID: '//trim(vname))
149 call nf90_err(nf90_get_var(ncid, varid, a2d), 'get variable: '//trim(vname))
150 call nf90_err(nf90_close(ncid), 'close: '//trim(gridfile))
151 kmt(:) = reshape(a2d, (/nxt*nyt/))
152 deallocate(a2d)
153 end if
154
155 ! -----------------------------------------------------------------------------
156 ! get the 3rd (vertical or ncat) dimension and variable attributes for the
157 ! ocean file
158 ! -----------------------------------------------------------------------------
159
160 call nf90_err(nf90_open(trim(input_file), nf90_nowrite, ncid), &
161 'open: '//trim(input_file))
162 if (do_ocnprep) then
163 call nf90_err(nf90_inq_dimid(ncid, 'Layer', varid), &
164 'get dimension Id: Layer'//trim(input_file))
165 call nf90_err(nf90_inquire_dimension(ncid, varid, len=nlevs), &
166 'get dimension Id: Layer'//trim(input_file))
167 else
168 call nf90_err(nf90_inq_dimid(ncid, 'ncat', varid), &
169 'get dimension Id: ncat'//trim(input_file))
170 call nf90_err(nf90_inquire_dimension(ncid, varid, len=nlevs), &
171 'get dimension Id: ncat'//trim(input_file))
172 endif
173 do n = 1,nvalid
174 if (do_ocnprep) then
175 if (trim(outvars(n)%var_name) .eq. 'eta')then
176 outvars(n)%long_name = 'Interface height'
177 outvars(n)%units = 'm'
178 else
179 call nf90_err(nf90_inq_varid(ncid, trim(outvars(n)%var_name), varid), &
180 'get variable Id: '//trim(outvars(n)%var_name))
181 call nf90_err(nf90_get_att(ncid, varid, 'long_name', outvars(n)%long_name), &
182 'get variable attribute: long_name '//trim(outvars(n)%var_name))
183 call nf90_err(nf90_get_att(ncid, varid, 'units', outvars(n)%units), &
184 'get variable attribute: units '//trim(outvars(n)%var_name) )
185 end if
186 end if
187 end do
188 call nf90_err(nf90_close(ncid), 'close: '//trim(input_file))
189
190 if (debug) then
191 do n = 1,nvalid
192 write(logunit,'(i4,a14,i4,a10,3(a6),a2)')n,' '//trim(outvars(n)%var_name)// &
193 ', ', outvars(n)%var_dimen,', '//trim(outvars(n)%var_remapmethod), &
194 ', '//trim(outvars(n)%var_grid), ', '//trim(outvars(n)%var_pair), &
195 ', '//trim(outvars(n)%var_pair_grid)
196 end do
197 end if
198
199 ! -----------------------------------------------------------------------------
200 ! get the masking variable for ocean 3-d remapping and create the mask
201 ! -----------------------------------------------------------------------------
202
203 if (do_ocnprep) then
204 allocate(eta(nlevs,nxt*nyt)); eta=0.0
205 call calc_eta(trim(input_file),(/nxt,nyt,nlevs/),bathysrc)
206
207 allocate(mask3d(nlevs,nxt*nyt)); mask3d = 0.0
208 call getfield(trim(input_file), trim(maskvar), dims=(/nxt,nyt,nlevs/), field=mask3d)
209 where(mask3d .le. real(hmin,4))mask3d = hmin
210
211 where(mask3d .le. hmin)mask3d = maskspval
212 where(mask3d .ne. maskspval)mask3d = 1.0
213
214 if (debug) then
215 call dumpnc(trim(ftype)//'.'//trim(fsrc)//'.eta.nc', 'eta', &
216 dims=(/nxt,nyt,nlevs/), field=eta)
217 call dumpnc(trim(ftype)//'.'//trim(fsrc)//'.mask3d.nc', 'mask3d', &
218 dims=(/nxt,nyt,nlevs/), field=mask3d)
219 end if
220 end if
221
222 ! -----------------------------------------------------------------------------
223 ! create packed arrays for mapping and remap arrays to the destination grid
224 ! -----------------------------------------------------------------------------
225
226 call setup_packing(nvalid,outvars)
227
228 ! 2D bilin
229 if (allocated(bilin2d)) then
230 call packarrays(trim(input_file), trim(wgtsdir)//fsrc(3:5)//'/', &
231 cos(angsrc), sin(angsrc), b2d, dims=(/nxt,nyt/), nflds=nbilin2d, fields=bilin2d)
232 rgb2d = 0.0
233 call remaprh(src_field=bilin2d, dst_field=rgb2d,rc=rc)
234 if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
235
236 if (debug) then
237 write(logunit,'(a)')'remap 2D fields bilinear with RH '
238 write(logunit,'(a)')'packed min/max values, mapped min/max values'
239 do n = 1,nbilin2d
240 write(logunit,'(i4,a14,3(a2,a6),4g14.4)')n,' '// &
241 trim(b2d(n)%var_name),' ',trim(b2d(n)%var_grid),' ', &
242 trim(b2d(n)%var_pair),' ',trim(b2d(n)%var_pair_grid), &
243 minval(bilin2d(n,:)), maxval(bilin2d(n,:)), &
244 minval(rgb2d(n,:)), maxval(rgb2d(n,:))
245 end do
246 call dumpnc(trim(ftype)//'.'//trim(fsrc)//'.bilin2d.nc', 'bilin2d', &
247 dims=(/nxt,nyt/), nflds=nbilin2d, field=bilin2d)
248 call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgbilin2d.nc', 'rgbilin2d', &
249 dims=(/nxr,nyr/), nflds=nbilin2d, field=rgb2d)
250 end if
251 end if
252
253 ! 2D conserv
254 if (allocated(consd2d)) then
255 call packarrays(trim(input_file), trim(wgtsdir)//fsrc(3:5)//'/', &
256 cos(angsrc), sin(angsrc), c2d, dims=(/nxt,nyt/), nflds=nconsd2d, fields=consd2d)
257 rgc2d = 0.0
258 call remaprh(src_field=consd2d, dst_field=rgc2d,rc=rc)
259 if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
260
261 if (debug) then
262 write(logunit,'(a)')'remap 2D fields conserv with RH '
263 write(logunit,'(a)')'packed min/max values, mapped min/max values'
264 do n = 1,nconsd2d
265 write(logunit,'(i4,a14,3(a2,a6),4g14.4)')n,' '// &
266 trim(c2d(n)%var_name),' ', trim(c2d(n)%var_grid),' ', &
267 trim(c2d(n)%var_pair),' ', trim(c2d(n)%var_pair_grid), &
268 minval(consd2d(n,:)), maxval(consd2d(n,:)), &
269 minval(rgc2d(n,:)), maxval(rgc2d(n,:))
270 end do
271 call dumpnc(trim(ftype)//'.'//trim(fsrc)//'.consd2d.nc', 'consd2d', &
272 dims=(/nxt,nyt/), nflds=nconsd2d, field=consd2d)
273 call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgconsd2d.nc', 'rgconsd2d', &
274 dims=(/nxr,nyr/), nflds=nconsd2d, field=rgc2d)
275 end if
276 end if
277
278 ! 3D bilin
279 if (allocated(bilin3d))then
280 call packarrays(trim(input_file), trim(wgtsdir)//fsrc(3:5)//'/', &
281 cos(angsrc), sin(angsrc), b3d, dims=(/nxt,nyt,nlevs/), nflds=nbilin3d, fields=bilin3d)
282
283 if (.not. do_ocnprep) then
284 ! -----------------------------------------------------------------------------
285 ! QC the source ice files after packing. Not every possible inconsistency is
286 ! checked but these are those known to create issues in the coupled model
287 ! -----------------------------------------------------------------------------
288
289 do n = 1,nbilin3d
290 if (trim(b3d(n)%var_name) == 'aicen')idx1 = n
291 if (trim(b3d(n)%var_name) == 'vicen')idx2 = n
292 if (trim(b3d(n)%var_name) == 'vsnon')idx3 = n
293 if (trim(b3d(n)%var_name) == 'Tsfcn')idx4 = n
294 end do
295
296 ! remove land values, if they exist
297 call zero_out_land_ice(kmt, bilin3d(idx1,:,:), icnt)
298 write(logunit, '(a,i8,a)')'removed ',icnt,' locations of '//trim(b3d(idx1)%var_name)//' from land'
299 call zero_out_land_ice(kmt, bilin3d(idx2,:,:), icnt)
300 write(logunit, '(a,i8,a)')'removed ',icnt,' locations of '//trim(b3d(idx2)%var_name)//' from land'
301 call zero_out_land_ice(kmt, bilin3d(idx3,:,:), icnt)
302 write(logunit, '(a,i8,a)')'removed ',icnt,' locations of '//trim(b3d(idx3)%var_name)//' from land'
303 call zero_out_land_ice(kmt, bilin3d(idx4,:,:), icnt)
304 write(logunit, '(a,i8,a)')'removed ',icnt,' locations of '//trim(b3d(idx4)%var_name)//' from land'
305
306 ! remove phantom-ice (aice=0, vice or vsno /= 0)
307 call zero_out_phantom_ice(bilin3d(idx1,:,:), bilin3d(idx2,:,:), icnt)
308 write(logunit, '(a,i8,a)')'removed ',icnt,' locations of phantom '//trim(b3d(idx2)%var_name)
309 call zero_out_phantom_ice(bilin3d(idx1,:,:), bilin3d(idx3,:,:), icnt)
310 write(logunit, '(a,i8,a)')'removed ',icnt,' locations of phantom '//trim(b3d(idx3)%var_name)
311 ! remove phantom-ice (vicen=0, aicen /=0); do not QC vsnon=0, aicen /=0)
312 call zero_out_phantom_ice(bilin3d(idx2,:,:), bilin3d(idx1,:,:), icnt)
313 write(logunit, '(a,i8,a)')'removed ',icnt,' locations of phantom '//trim(b3d(idx1)%var_name)
314 deallocate(kmt)
315 end if
316
317 rgb3d = 0.0
318 do k = 1,nlevs
319 if (do_ocnprep) then
320 call remaprh(k,src_field=bilin3d(:,k,:), dst_field=rgb3d(:,k,:), hmask=mask3d(k,:),rc=rc)
321 if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
322 else
323 call remaprh(src_field=bilin3d(:,k,:), dst_field=rgb3d(:,k,:),rc=rc)
324 if (chkerr(rc,__line__,u_file_u)) call esmf_finalize(endflag=esmf_end_abort)
325 end if
326 end do
327 if (do_ocnprep) then
328 call vfill()
329 end if
330
331 if (debug) then
332 write(logunit,'(a)')'remap 3D fields bilinear with RH '
333 write(logunit,'(a)')'packed min/max values,mapped min/max values'
334 do n = 1,nbilin3d
335 write(logunit,'(i4,a14,3(a2,a6),4g14.4)')n,' '// &
336 trim(b3d(n)%var_name),' ', trim(b3d(n)%var_grid),' ', &
337 trim(b3d(n)%var_pair),' ', trim(b3d(n)%var_pair_grid), &
338 minval(bilin3d(n,:,:)), maxval(bilin3d(n,:,:)), &
339 minval(rgb3d(n,:,:)), maxval(rgb3d(n,:,:))
340 end do
341 call dumpnc(trim(ftype)//'.'//trim(fsrc)//'.bilin3d.nc', 'bilin3d', &
342 dims=(/nxt,nyt,nlevs/), nk=nlevs, nflds=nbilin3d, field=bilin3d)
343 call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgbilin3d.nc', 'rgbilin3d', &
344 dims=(/nxr,nyr,nlevs/), nk=nlevs, nflds=nbilin3d, field=rgb3d)
345 end if
346 end if
347
348 !------------------------------------------------------------------------------
349 ! rotate on Ct from EN->IJ and remap back to native staggers
350 !------------------------------------------------------------------------------
351
352 if (allocated(bilin2d)) then
353 call rotremap(trim(wgtsdir)//fdst(3:5)//'/', b2d, cos(angdst), sin(angdst), &
354 dims=(/nxr,nyr/), nflds=nbilin2d, fields=rgb2d)
355 if (debug) then
356 call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgbilin2d.ij.nc', 'rgbilin2d', &
357 dims=(/nxr,nyr/), nflds=nbilin2d, field=rgb2d)
358 end if
359 end if
360 if (allocated(consd2d)) then
361 call rotremap(trim(wgtsdir)//fdst(3:5)//'/', c2d, cos(angdst), sin(angdst), &
362 dims=(/nxr,nyr/), nflds=nconsd2d, fields=rgc2d)
363 if (debug) then
364 call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgbilin2d.ij.nc', 'rgbilin2d', &
365 dims=(/nxr,nyr/), nflds=nconsd2d, field=rgc2d)
366 end if
367 end if
368 if (allocated(bilin3d)) then
369 call rotremap(trim(wgtsdir)//fdst(3:5)//'/', b3d, cos(angdst), sin(angdst), &
370 dims=(/nxr,nyr,nlevs/), nflds=nbilin3d, fields=rgb3d)
371 if (debug) then
372 call dumpnc(trim(ftype)//'.'//trim(fdst)//'.rgbilin3d.ij.nc', 'rgbilin3d', &
373 dims=(/nxr,nyr,nlevs/), nk=nlevs, nflds=nbilin3d, field=rgb3d)
374 end if
375 end if
376
377 ! -----------------------------------------------------------------------------
378 ! write the mapped fields
379 ! -----------------------------------------------------------------------------
380
381 allocate(out2d(nxr,nyr)); out2d = 0.0
382 allocate(out3d(nxr,nyr,nlevs)); out3d = 0.0
383
384 fout = trim(ftype)//'.'//trim(fdst)//'.nc'
385 if (debug) write(logunit, '(a)')'output file: '//trim(fout)
386
387 if (do_ocnprep) then
388 call setup_ocnrestart(trim(input_file),trim(fout),bathydst)
389 else
390 call setup_icerestart(trim(input_file),trim(fout))
391 end if
392
393 call nf90_err(nf90_open(trim(fout), nf90_write, ncid), 'write: '//trim(fout))
394 if (allocated(rgb2d)) then
395 do n = 1,nbilin2d
396 out2d(:,:) = reshape(rgb2d(n,:), (/nxr,nyr/))
397 if (b2d(n)%var_grid(1:2) == 'Bu') out2d(:,nyr) = out2d(:,nyr-1)
398 vname = trim(b2d(n)%var_name)
399 call nf90_err(nf90_inq_varid(ncid, vname, varid), 'get variable Id: '//vname)
400 call nf90_err(nf90_put_var(ncid, varid, out2d), 'put variable: '//vname)
401 end do
402 end if
403 if (allocated(rgc2d)) then
404 do n = 1,nconsd2d
405 out2d(:,:) = reshape(rgc2d(n,:), (/nxr,nyr/))
406 vname = trim(c2d(n)%var_name)
407 call nf90_err(nf90_inq_varid(ncid, vname, varid), 'get variable Id: '//vname)
408 call nf90_err(nf90_put_var(ncid, varid, out2d), 'put variable: '//vname)
409 end do
410 end if
411 if (allocated(rgb3d)) then
412 do n = 1,nbilin3d
413 do k = 1,nlevs
414 out3d(:,:,k) = reshape(rgb3d(n,k,:), (/nxr,nyr/))
415 end do
416 if (b3d(n)%var_grid(1:2) == 'Cv') out3d(:,nyr,:) = out3d(:,nyr-1,:)
417 vname = trim(b3d(n)%var_name)
418 call nf90_err(nf90_inq_varid(ncid, vname, varid), 'get variable Id: '//vname)
419 call nf90_err(nf90_put_var(ncid, varid, out3d), 'put variable: '//vname)
420 end do
421 end if
422 call nf90_err(nf90_close(ncid), 'close: '// trim(fout))
423 write(logunit,'(a)')trim(fout)//' done'
424 stop
425
426end program ocniceprep
utils_esmf_mod::remaprh
Definition utils_esmf_mod.F90:26
utils_esmf_mod::rotremap
Definition utils_esmf_mod.F90:33
utils_mod::dumpnc
Definition utils_mod.F90:38
utils_mod::getfield
Definition utils_mod.F90:17
utils_mod::packarrays
Definition utils_mod.F90:22
utils_mod::remap
Definition utils_mod.F90:27
ocniceprep
program ocniceprep
Read either a MOM6 or CICE6 restart file at 1/4deg tripole resolution and remap the required fields t...
Definition ocniceprep.F90:18