22 use gengrid_kinds,
only: cl, cs, dbl_kind, real_kind, int_kind
23 use angles,
only: find_angq, find_ang
24 use vertices,
only: fill_vertices, fill_bottom, fill_top
25 use mapped_mask,
only: make_frac_land
26 use postwgts,
only: make_postwgts
27 use tripolegrid,
only: write_tripolegrid
28 use cicegrid,
only: write_cicegrid
29 use scripgrid,
only: write_scripgrid
30 use topoedits,
only: add_topoedits, apply_topoedits
31 use charstrings,
only: logmsg, res, dirsrc, dirout, atmres, fv3dir, editsfile
32 use charstrings,
only: maskfile, maskname, topofile, toponame, editsfile, staggerlocs, cdate, history
33 use debugprint,
only: checkseam, checkxlatlon, checkpoint
39 real(dbl_kind) :: dxt, dyt
41 real(kind=dbl_kind),
parameter :: pi = 3.14159265358979323846_dbl_kind
42 real(kind=dbl_kind),
parameter :: deg2rad = pi/180.0_dbl_kind
44 real(real_kind),
allocatable,
dimension(:,:) :: ww3dpth
45 integer(int_kind),
allocatable,
dimension(:,:) :: ww3mask
47 character(len=CL) :: fsrc, fdst, fwgt
48 character(len= 2) :: cstagger
50 integer :: rc,ncid,id,xtype
51 integer :: i,j,k,i2,j2
53 integer :: localpet, npet
54 logical :: fexist = .false.
56 type(esmf_regridmethod_flag) :: method
60 character(len= 6) :: i4fmt =
'(i4.4)' 61 character(len=CS) :: form1
62 character(len=CS) :: form2
63 character(len= 6) :: cnx
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)
75 if (localpet == 0) mastertask=.true.
77 print *,npet,
' More than one task specified; Aborting ' 78 call esmf_finalize(endflag=esmf_end_abort)
85 call read_inputnml(
'grid.nml')
88 print
'(a,2i6)',
' output grid requested ',ni,nj
89 print
'(a,2i6)',
' supergrid size used ', nx,ny
90 print
'(a)',
' output grid tag '//trim(res)
91 print
'(a)',
' supergrid source directory '//trim(dirsrc)
92 print
'(a)',
' output grid directory '//trim(dirout)
93 print
'(a)',
' atm resolution '//trim(atmres)
94 print
'(a,i6)',
' fv3 tile grid size ',npx
95 print
'(a)',
' atm mosaic directory '//trim(fv3dir)
96 print
'(a)',
' MOM6 topography file '//trim(topofile)
97 print
'(a)',
' MOM6 edits file '//trim(editsfile)
98 print *,
'editmask flag ',editmask
99 print *,
'debug flag ',debug
100 print *,
'do_postwgts flag ',do_postwgts
106 call esmf_logwrite(
"Starting gen_fixgrid", esmf_logmsg_info)
111 ivertcu = ivertct + 1; jvertcu = jvertct + 0
112 ivertcv = ivertct + 0; jvertcv = jvertct + 1
113 ivertbu = ivertct + 1; jvertbu = jvertct + 1
116 print
'(a8,4i6)',
'iVertCt ',(ivertct(i),i=1,4)
117 print
'(a8,4i6)',
'jVertCt ',(jvertct(i),i=1,4)
119 print
'(a8,4i6)',
'iVertCu ',(ivertcu(i),i=1,4)
120 print
'(a8,4i6)',
'jVertCu ',(jvertcu(i),i=1,4)
122 print
'(a8,4i6)',
'iVertCv ',(ivertcv(i),i=1,4)
123 print
'(a8,4i6)',
'jVertCv ',(jvertcv(i),i=1,4)
125 print
'(a8,4i6)',
'iVertBu ',(ivertbu(i),i=1,4)
126 print
'(a8,4i6)',
'jVertBu ',(jvertbu(i),i=1,4)
130 latct_vert = -9999.0 ; lonct_vert = -9999.0
131 latcu_vert = -9999.0 ; loncu_vert = -9999.0
132 latcv_vert = -9999.0 ; loncv_vert = -9999.0
133 latbu_vert = -9999.0 ; lonbu_vert = -9999.0
139 fsrc = trim(dirsrc)//
'/'//trim(maskfile)
141 rc = nf90_open(fsrc, nf90_nowrite, ncid)
143 print
'(a)',
'reading ocean mask from '//trim(fsrc)
144 if(rc .ne. 0)print
'(a)',
'nf90_open = '//trim(nf90_strerror(rc))
147 wet4 = 0.0; wet8 = 0.0
148 rc = nf90_inq_varid(ncid, trim(maskname), id)
149 rc = nf90_inquire_variable(ncid, id, xtype=xtype)
150 if(xtype .eq. 5)rc = nf90_get_var(ncid, id, wet4)
151 if(xtype .eq. 6)rc = nf90_get_var(ncid, id, wet8)
152 rc = nf90_close(ncid)
154 if(xtype.eq. 6)wet4 =
real(wet8,4)
163 fsrc = trim(dirsrc)//
'/'//trim(topofile)
165 rc = nf90_open(fsrc, nf90_nowrite, ncid)
167 print
'(a)',
'reading ocean topography from '//trim(fsrc)
168 if(rc .ne. 0)print
'(a)',
'nf90_open = '//trim(nf90_strerror(rc))
172 rc = nf90_inq_varid(ncid, trim(toponame), id)
173 rc = nf90_inquire_variable(ncid, id, xtype=xtype)
174 if(xtype .eq. 5)rc = nf90_get_var(ncid, id, dp4)
175 if(xtype .eq. 6)rc = nf90_get_var(ncid, id, dp8)
176 rc = nf90_close(ncid)
178 if(xtype.eq. 6)dp4 =
real(dp8,4)
188 if(trim(editsfile) ==
'none')
then 189 print
'(a)',
'Need a valid editsfile to make mask edits ' 193 fsrc = trim(dirsrc)//
'/'//trim(editsfile)
194 fdst = trim(dirout)//
'/'//
'ufs.'//trim(editsfile)
195 call add_topoedits(fsrc,fdst)
204 fsrc = trim(dirsrc)//
'/'//trim(editsfile)
205 if(editmask)fsrc = trim(dirout)//
'/'//
'ufs.'//trim(editsfile)
206 call apply_topoedits(fsrc)
212 fsrc = trim(dirsrc)//
'/'//
'ocean_hgrid.nc' 214 rc = nf90_open(fsrc, nf90_nowrite, ncid)
216 print
'(a)',
'reading supergrid from '//trim(fsrc)
217 if(rc .ne. 0)print
'(a)',
'nf90_open = '//trim(nf90_strerror(rc))
220 rc = nf90_inq_varid(ncid,
'x', id)
221 rc = nf90_get_var(ncid, id, x)
223 rc = nf90_inq_varid(ncid,
'y', id)
224 rc = nf90_get_var(ncid, id, y)
226 rc = nf90_inq_varid(ncid,
'dx', id)
227 rc = nf90_get_var(ncid, id, dx)
229 rc = nf90_inq_varid(ncid,
'dy', id)
230 rc = nf90_get_var(ncid, id, dy)
232 rc = nf90_close(ncid)
235 sg_maxlat = maxval(y)
251 ulon(i,j) = x(i2,j2)*deg2rad
252 ulat(i,j) = y(i2,j2)*deg2rad
254 angle(i,j) = -angq(i2,j2)
256 htn(i,j) = (dx(i2-1,j2) + dx(i2,j2))*100._dbl_kind
257 hte(i,j) = (dy(i2,j2-1) + dy(i2,j2))*100._dbl_kind
259 lonbu(i,j) = x(i2,j2)
260 latbu(i,j) = y(i2,j2)
262 lonct(i,j) = x(i2-1,j2-1)
263 loncu(i,j) = x(i2, j2-1)
264 loncv(i,j) = x(i2-1,j2 )
266 latct(i,j) = y(i2-1,j2-1)
267 latcu(i,j) = y(i2, j2-1)
268 latcv(i,j) = y(i2-1,j2 )
270 dxt = dx(i2-1,j2-1) + dx(i2,j2-1)
271 dyt = dy(i2-1,j2-1) + dy(i2-1,j2)
272 areact(i,j) = dxt*dyt
283 print *,
'ANGLET ',minval(anglet),maxval(anglet)
284 print *,
'ANGLE ',minval(angle),maxval(angle)
296 write(logmsg,
'(a,2e12.5)')
'min vals of hte at folds ', &
297 minval(hte(ni/2,:)),minval(hte(ni,:))
298 print
'(a)',trim(logmsg)
302 if(hte(ii,j) .le. 1.0)hte(ii,j) = 0.5*(hte(ii-1,j) + hte(ii+1,j))
304 if(hte(ii,j) .le. 1.0)hte(ii,j) = 0.5*(hte(ii-1,j) + hte( 1,j))
307 write(logmsg,
'(a,2e12.5)')
'min vals of hte at folds ', &
308 minval(hte(ni/2,:)),minval(hte(ni,:))
309 print
'(a)',trim(logmsg)
316 where(lonct .lt. 0.0)lonct = lonct + 360._dbl_kind
317 where(loncu .lt. 0.0)loncu = loncu + 360._dbl_kind
318 where(loncv .lt. 0.0)loncv = loncv + 360._dbl_kind
319 where(lonbu .lt. 0.0)lonbu = lonbu + 360._dbl_kind
330 if(latbu(i,j) .eq. sg_maxlat)ipole(1) = i
333 if(latbu(i,j) .eq. sg_maxlat)ipole(2) = i
336 write(logmsg,
'(a,2i6,2f12.2)')
'poles found at i = ',ipole,latbu(ipole(1),nj), &
338 print
'(a)',trim(logmsg)
341 if(mastertask .and. debug)
call checkseam
344 i2 = ipole(2)+(ipole(1)-i)+1
345 xlonct(i) = lonct(i2,nj)
346 xlatct(i) = latct(i2,nj)
350 i2 = ipole(2)+(ipole(1)-i)
352 xloncu(i) = loncu(i2,nj)
353 xlatcu(i) = latcu(i2,nj)
356 if(mastertask .and. debug)
call checkxlatlon
360 dlatbu(i) = latbu(i,1) + 2.0*(latcu(i,1) - latbu(i,1))
361 dlatcv(i) = latct(i,1) + 2.0*(latct(i,1) - latcv(i,1))
369 call fill_vertices(2,nj , ivertct,jvertct, latbu,lonbu, latct_vert,lonct_vert)
370 call fill_bottom(ivertct,jvertct, latbu,lonbu, latct_vert,lonct_vert,dlatbu)
372 call fill_vertices(2,nj , ivertcu,jvertcu, latcv,loncv, latcu_vert,loncu_vert)
373 call fill_bottom(ivertcu,jvertcu, latcv,loncv, latcu_vert,loncu_vert,dlatcv)
376 call fill_vertices(1,nj-1, ivertcv,jvertcv, latcu,loncu, latcv_vert,loncv_vert)
377 call fill_top(ivertcv,jvertcv, latcu,loncu, latcv_vert,loncv_vert, xlatcu, xloncu)
379 call fill_vertices(1,nj-1, ivertbu,jvertbu, latct,lonct, latbu_vert,lonbu_vert)
380 call fill_top(ivertbu,jvertbu, latct,lonct, latbu_vert,lonbu_vert, xlatct, xlonct)
382 if(mastertask .and. debug)
call checkpoint
384 if(minval(latct_vert) .lt. -1.e3)stop
385 if(minval(lonct_vert) .lt. -1.e3)stop
386 if(minval(latcu_vert) .lt. -1.e3)stop
387 if(minval(loncu_vert) .lt. -1.e3)stop
388 if(minval(latcv_vert) .lt. -1.e3)stop
389 if(minval(loncv_vert) .lt. -1.e3)stop
390 if(minval(latbu_vert) .lt. -1.e3)stop
391 if(minval(lonbu_vert) .lt. -1.e3)stop
392 deallocate(xlonct, xlatct, xloncu, xlatcu, dlatbu, dlatcv)
399 call date_and_time(date=cdate)
400 history =
'created on '//trim(cdate)//
' from '//trim(fsrc)
403 fdst = trim(dirout)//
'/'//
'tripole.mx'//trim(res)//
'.nc' 404 call write_tripolegrid(trim(fdst))
407 fdst = trim(dirout)//
'/'//
'grid_cice_NEMS_mx'//trim(res)//
'.nc' 408 call write_cicegrid(trim(fdst))
409 deallocate(ulon, ulat, htn, hte)
413 cstagger = trim(staggerlocs(k))
414 fdst = trim(dirout)//
'/'//trim(cstagger)//
'.mx'//trim(res)//
'_SCRIP.nc' 416 logmsg =
'creating SCRIP file '//trim(fdst)
417 print
'(a)',trim(logmsg)
419 call write_scripgrid(trim(fdst),trim(cstagger))
421 deallocate(latcv_vert, loncv_vert)
422 deallocate(latcu_vert, loncu_vert)
423 deallocate(latbu_vert, lonbu_vert)
427 cstagger = trim(staggerlocs(1))
428 fdst= trim(dirout)//
'/'//trim(cstagger)//
'.mx'//trim(res)//
'_SCRIP_land.nc' 430 logmsg =
'creating SCRIP file '//trim(fdst)
431 print
'(a)',trim(logmsg)
433 call write_scripgrid(trim(fdst),trim(cstagger),imask=int(wet4))
434 deallocate(latct_vert, lonct_vert)
442 write(form1,
'(a)')
'('//trim(cnx)//
'f14.8)' 443 write(form2,
'(a)')
'('//trim(cnx)//
'i2)' 445 allocate(ww3mask(1:ni,1:nj)); ww3mask = wet4
446 allocate(ww3dpth(1:ni,1:nj)); ww3dpth = dp4
448 where(latct .ge. maximum_lat)ww3mask = 3
452 open(unit=21,file=trim(dirout)//
'/'//
'ww3.mx'//trim(res)//
'_x.inp',form=
'formatted')
453 open(unit=22,file=trim(dirout)//
'/'//
'ww3.mx'//trim(res)//
'_y.inp',form=
'formatted')
454 open(unit=23,file=trim(dirout)//
'/'//
'ww3.mx'//trim(res)//
'_bottom.inp',form=
'formatted')
455 open(unit=24,file=trim(dirout)//
'/'//
'ww3.mx'//trim(res)//
'_mapsta.inp',form=
'formatted')
457 open(unit=25,file=trim(dirout)//
'/'//
'ww3.mx'//trim(res)//
'_obstr.inp',form=
'formatted')
460 write( 21,trim(form1))lonct(:,j)
461 write( 22,trim(form1))latct(:,j)
464 write( 23,trim(form1))ww3dpth(:,j)
465 write( 24,trim(form2))ww3mask(:,j)
467 write( 25,trim(form2))ww3mask(:,j)*0
468 write( 25,trim(form2))ww3mask(:,j)*0
471 close(21);
close(22);
close(23);
close(24);
close(25)
472 deallocate(ww3mask);
deallocate(ww3dpth)
473 deallocate(wet4, wet8)
481 method=esmf_regridmethod_conserve
482 fsrc = trim(dirout)//
'/'//
'Ct.mx'//trim(res)//
'_SCRIP_land.nc' 483 fdst = trim(fv3dir)//
'/'//trim(atmres)//
'/'//trim(atmres)//
'_mosaic.nc' 484 fwgt = trim(dirout)//
'/'//
'Ct.mx'//trim(res)//
'.to.'//trim(atmres)//
'.nc' 486 logmsg =
'creating weight file '//trim(fwgt)
487 print
'(a)',trim(logmsg)
490 call esmf_regridweightgen(srcfile=trim(fsrc),dstfile=trim(fdst), &
491 weightfile=trim(fwgt), regridmethod=method, &
492 unmappedaction=esmf_unmappedaction_ignore, &
493 ignoredegenerate=.true., &
494 tilefilepath=trim(fv3dir)//
'/'//trim(atmres)//
'/', rc=rc)
495 if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
496 line=__line__, file=__file__))
call esmf_finalize(endflag=esmf_end_abort)
499 logmsg =
'creating mapped ocean mask for '//trim(atmres)
500 print
'(a)',trim(logmsg)
502 call make_frac_land(trim(fsrc), trim(fwgt))
510 if(trim(res) .ne.
'025')
then 511 fsrc = trim(dirout)//
'/'//
'Ct.mx025_SCRIP.nc' 512 inquire(file=trim(fsrc), exist=fexist)
514 method=esmf_regridmethod_nearest_stod
515 fdst = trim(dirout)//
'/'//
'Ct.mx'//trim(res)//
'_SCRIP.nc' 516 fwgt = trim(dirout)//
'/'//
'tripole.mx025.Ct.to.mx'//trim(res)//
'.Ct.neareststod.nc' 518 logmsg =
'creating weight file '//trim(fwgt)
519 print
'(a)',trim(logmsg)
522 call esmf_regridweightgen(srcfile=trim(fsrc),dstfile=trim(fdst), &
523 weightfile=trim(fwgt), regridmethod=method, &
524 ignoredegenerate=.true., &
525 unmappedaction=esmf_unmappedaction_ignore, rc=rc)
526 if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
527 line=__line__, file=__file__))
call esmf_finalize(endflag=esmf_end_abort)
530 logmsg =
'ERROR: '//trim(fsrc)//
' is required to generate tripole:triple weights' 531 print
'(a)',trim(logmsg)
541 if(do_postwgts)
call make_postwgts
547 deallocate(x,y, angq, dx, dy, xsgp1, ysgp1)
548 deallocate(areact, anglet, angle)
549 deallocate(latct, lonct)
550 deallocate(latcv, loncv)
551 deallocate(latcu, loncu)
552 deallocate(latbu, lonbu)
program gen_fixgrid
Generate fixed grid files required for coupled model using the MOM6 super grid file and ocean mask fi...