gen_fixgrid.F90

Go to the documentation of this file.

 
 
program gen_fixgrid
 
  use esmf
 
  use grdvars
  use inputnml
  use gengrid_kinds,     only: cl, cs, dbl_kind, real_kind, int_kind
  use angles,            only: find_ang, find_angq, find_angchk
  use vertices,          only: fill_vertices, fill_bottom, fill_top
  use mapped_mask,       only: make_frac_land
  use postwgts,          only: make_postwgts
  use tripolegrid,       only: write_tripolegrid
  use cicegrid,          only: write_cicegrid
  use scripgrid,         only: write_scripgrid
  use topoedits,         only: add_topoedits, apply_topoedits
  use charstrings,       only: logmsg, res, atmres, dirsrc, dirout, fv3dir, editsfile
  use charstrings,       only: maskfile, maskname, topofile, toponame, editsfile, staggerlocs, cdate, history
  use debugprint,        only: checkseam, checkxlatlon, checkpoint
  use netcdf
 
  implicit none
 
  ! local variables
  real(dbl_kind) :: dxt, dyt
 
  real(kind=dbl_kind), parameter :: pi = 3.14159265358979323846_dbl_kind
  real(kind=dbl_kind), parameter :: deg2rad = pi/180.0_dbl_kind
 
  real(real_kind),   allocatable, dimension(:,:) :: ww3dpth
  integer(int_kind), allocatable, dimension(:,:) :: ww3mask
 
  character(len=CL) :: fsrc, fdst, fwgt
  character(len= 2) :: cstagger
 
  integer :: rc,ncid,id,xtype
  integer :: i,j,k,n,i2,j2
  integer :: ii
  integer :: localpet, npet
  logical :: fexist = .false.
 
  type(esmf_regridmethod_flag) :: method
  type(esmf_vm) :: vm
 
  !WW3 file format for mod_def generation
  character(len= 6) :: i4fmt = '(i4.4)'
  character(len=CS) :: form1
  character(len=CS) :: form2
  character(len= 6) :: cnx
 
  !-------------------------------------------------------------------------
  ! Initialize esmf environment.
  !-------------------------------------------------------------------------
 
  call esmf_vmgetglobal(vm, rc=rc)
  call esmf_initialize(vm=vm, logkindflag=esmf_logkind_multi, rc=rc)
  call esmf_vmget(vm, localpet=localpet, pecount=npet, rc=rc)
  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
       line=__line__, file=__file__)) call esmf_finalize(endflag=esmf_end_abort)
  roottask = .false.
  if (localpet == 0) roottask =.true.
  if (npet /= 1) then
     print *,npet,' More than one task specified; Aborting '
     call esmf_finalize(endflag=esmf_end_abort)
  end if
 
  !---------------------------------------------------------------------
  !
  !---------------------------------------------------------------------
 
  call read_inputnml('grid.nml')
 
  print '(a,2i6)',' output grid requested ',ni,nj
  print '(a,2i6)',' supergrid size used ', nx,ny
  print '(a)',' output grid tag '//trim(res)
  print '(a)',' supergrid source directory '//trim(dirsrc)
  print '(a)',' output grid directory '//trim(dirout)
  print '(a)',' atm mosaic directory '//trim(fv3dir)
  print '(a)',' MOM6 topography file '//trim(topofile)
  print '(a)',' MOM6 edits file '//trim(editsfile)
  print *,'editmask flag ',editmask
  print *,'debug flag ',debug
  print *,'do_postwgts flag ',do_postwgts
  print *
 
  call allocate_all
 
  call esmf_logwrite("Starting gen_fixgrid", esmf_logmsg_info)
  !---------------------------------------------------------------------
  ! set up the arrays to retrieve the vertices
  !---------------------------------------------------------------------
 
  ivertcu = ivertct + 1; jvertcu = jvertct + 0
  ivertcv = ivertct + 0; jvertcv = jvertct + 1
  ivertbu = ivertct + 1; jvertbu = jvertct + 1
 
  print '(a8,4i6)','iVertCt ',(ivertct(i),i=1,4)
  print '(a8,4i6)','jVertCt ',(jvertct(i),i=1,4)
  print *
  print '(a8,4i6)','iVertCu ',(ivertcu(i),i=1,4)
  print '(a8,4i6)','jVertCu ',(jvertcu(i),i=1,4)
  print *
  print '(a8,4i6)','iVertCv ',(ivertcv(i),i=1,4)
  print '(a8,4i6)','jVertCv ',(jvertcv(i),i=1,4)
  print *
  print '(a8,4i6)','iVertBu ',(ivertbu(i),i=1,4)
  print '(a8,4i6)','jVertBu ',(jvertbu(i),i=1,4)
  print *
 
  latct_vert = -9999.0 ; lonct_vert = -9999.0
  latcu_vert = -9999.0 ; loncu_vert = -9999.0
  latcv_vert = -9999.0 ; loncv_vert = -9999.0
  latbu_vert = -9999.0 ; lonbu_vert = -9999.0
 
  !---------------------------------------------------------------------
  ! read the MOM6 land mask
  !---------------------------------------------------------------------
 
  fsrc = trim(dirsrc)//'/'//trim(maskfile)
 
  rc = nf90_open(fsrc, nf90_nowrite, ncid)
  print '(a)', 'reading ocean mask from '//trim(fsrc)
  if(rc .ne. 0)print '(a)', 'nf90_open = '//trim(nf90_strerror(rc))
 
  wet4 = 0.0; wet8 = 0.0
  rc = nf90_inq_varid(ncid,  trim(maskname), id)
  rc = nf90_inquire_variable(ncid, id, xtype=xtype)
  if(xtype .eq. 5)rc = nf90_get_var(ncid,      id,  wet4)
  if(xtype .eq. 6)rc = nf90_get_var(ncid,      id,  wet8)
  rc = nf90_close(ncid)
 
  if(xtype.eq. 6)wet4 = real(wet8,4)
 
  !---------------------------------------------------------------------
  ! read the MOM6 depth file
  !---------------------------------------------------------------------
 
  fsrc = trim(dirsrc)//'/'//trim(topofile)
 
  rc = nf90_open(fsrc, nf90_nowrite, ncid)
  print '(a)', 'reading ocean topography from '//trim(fsrc)
  if(rc .ne. 0)print '(a)', 'nf90_open = '//trim(nf90_strerror(rc))
 
  dp4 = 0.0; dp8 = 0.0
  rc = nf90_inq_varid(ncid,  trim(toponame), id)
  rc = nf90_inquire_variable(ncid, id, xtype=xtype)
  if(xtype .eq. 5)rc = nf90_get_var(ncid,      id,  dp4)
  if(xtype .eq. 6)rc = nf90_get_var(ncid,      id,  dp8)
  rc = nf90_close(ncid)
 
  if(xtype.eq. 6)dp4 = real(dp8,4)
 
  if(editmask)then
     !---------------------------------------------------------------------
     ! apply topoedits run time mask changes if required for this config
     ! this will create a modified topoedits file which accounts for any
     ! land mask changes created at run time by MOM6
     !---------------------------------------------------------------------
 
     if(trim(editsfile)  == 'none')then
        print '(a)', 'Need a valid editsfile to make mask edits '
        call abort()
     end if
     inquire(file=trim(dirsrc)//'/'//trim(editsfile),exist=fexist)
     if (.not. fexist) then
        print '(a)', 'Required topoedits file '//trim(editsfile) &
             //'for land mask changes is missing '
        call abort()
     end if
 
     fsrc = trim(dirsrc)//'/'//trim(editsfile)
     fdst = trim(dirout)//'/'//'ufs.'//trim(editsfile)
     call add_topoedits(fsrc,fdst)
  endif
 
  !---------------------------------------------------------------------
  ! MOM6 reads the depth file, applies the topo edits and then adjusts
  ! depth using masking_depth and min/max depth. This call mimics
  ! MOM6 routines apply_topography_edits_from_file and limit_topography
  ! If the the topoedits file has been modified to account for MOM6 run
  ! time land mask changes (above), then the depth will be created using
  ! this modified topoedits file
  !---------------------------------------------------------------------
 
  fsrc = trim(dirsrc)//'/'//trim(editsfile)
  if(editmask)fsrc = trim(dirout)//'/'//'ufs.'//trim(editsfile)
 
  if (trim(editsfile) /= 'none') then
     inquire(file=trim(fsrc),exist=fexist)
     if (.not. fexist) then
        print '(a)', 'Required topoedits file '//trim(fsrc)//' is missing '
        call abort()
     end if
  end if
  call apply_topoedits(fsrc)
 
  !---------------------------------------------------------------------
  ! read MOM6 supergrid file
  !---------------------------------------------------------------------
 
  fsrc = trim(dirsrc)//'/'//'ocean_hgrid.nc'
 
  rc = nf90_open(fsrc, nf90_nowrite, ncid)
  print '(a)', 'reading supergrid from '//trim(fsrc)
  if(rc .ne. 0)print '(a)', 'nf90_open = '//trim(nf90_strerror(rc))
 
  rc = nf90_inq_varid(ncid, 'x', id)  !lon
  rc = nf90_get_var(ncid,    id,  x)
 
  rc = nf90_inq_varid(ncid, 'y', id)  !lat
  rc = nf90_get_var(ncid,    id,  y)
 
  rc = nf90_inq_varid(ncid, 'dx', id)
  rc = nf90_get_var(ncid,     id, dx)
 
  rc = nf90_inq_varid(ncid, 'dy', id)
  rc = nf90_get_var(ncid,     id, dy)
 
  rc = nf90_close(ncid)
  sg_maxlat = maxval(y)
  write(logmsg,'(a,f12.2)')'max lat in super grid ',maxval(y)
  print '(a)',trim(logmsg)
 
  !---------------------------------------------------------------------
  ! fill grid variables
  !---------------------------------------------------------------------
 
  do j = 1,nj
     do i = 1,ni
        i2 = 2*i ; j2 = 2*j
        !deg->rad
        ulon(i,j) =     x(i2,j2)*deg2rad
        ulat(i,j) =     y(i2,j2)*deg2rad
        !m->cm
        htn(i,j) = (dx(i2-1,j2) + dx(i2,j2))*100._dbl_kind
        hte(i,j) = (dy(i2,j2-1) + dy(i2,j2))*100._dbl_kind
        !deg
        lonbu(i,j) =     x(i2,j2)
        latbu(i,j) =     y(i2,j2)
        !deg
        lonct(i,j) =     x(i2-1,j2-1)
        loncu(i,j) =     x(i2,  j2-1)
        loncv(i,j) =     x(i2-1,j2  )
        !deg
        latct(i,j) =     y(i2-1,j2-1)
        latcu(i,j) =     y(i2,  j2-1)
        latcv(i,j) =     y(i2-1,j2  )
        !m2
        dxt = dx(i2-1,j2-1) + dx(i2,j2-1)
        dyt = dy(i2-1,j2-1) + dy(i2-1,j2)
        areact(i,j) = dxt*dyt
     enddo
  enddo
 
  !---------------------------------------------------------------------
  ! locate the ith index of the two poles on j=nj
  ! the corner points must lie on the pole
  !---------------------------------------------------------------------
 
  ipole = -1
  j = nj
  do i = 1,ni/2
     if(latbu(i,j) .eq. sg_maxlat)ipole(1) = i
  enddo
  do i = ni/2+1,ni
     if(latbu(i,j) .eq. sg_maxlat)ipole(2) = i
  enddo
  write(logmsg,'(a,2i6,2f12.2)')'poles found at i = ',ipole, latbu(ipole(1),nj), &
       latbu(ipole(2),nj)
  print '(a)',trim(logmsg)
 
  !---------------------------------------------------------------------
  ! find the angle on centers using the same procedure as MOM6
  !---------------------------------------------------------------------
 
  call find_ang((/1,ni/),(/1,nj/),lonbu,latbu,lonct,anglet)
  write(logmsg,'(a,2f12.2)')'ANGLET min,max: ',minval(anglet),maxval(anglet)
  print '(a)',trim(logmsg)
  write(logmsg,'(a,2f12.2)')'ANGLET edges i=1,i=ni: ',anglet(1,nj),anglet(ni,nj)
  print '(a)',trim(logmsg)
 
  xangct(:) = 0.0
  do i = 1,ni
     i2 = ipole(2)+(ipole(1)-i)+1
     xangct(i) = -anglet(i2,nj)       ! angle changes sign across seam
  end do
 
  !---------------------------------------------------------------------
  ! find the angle on corners using the same procedure as CICE6
  !---------------------------------------------------------------------
 
  call find_angq((/1,ni/),(/1,nj/),xangct,anglet,angle)
  angle(ni,:) = -angle(1,:)
  ! reverse angle for CICE
  angle = -angle
  write(logmsg,'(a,2f12.2)')'ANGLE min,max: ',minval(angle),maxval(angle)
  print '(a)',trim(logmsg)
  write(logmsg,'(a,2f12.2)')'ANGLE edges i=1,i=ni: ',angle(1,nj),angle(ni,nj)
  print '(a)',trim(logmsg)
 
  !---------------------------------------------------------------------
  ! check the Bu angle
  !---------------------------------------------------------------------
 
  call find_angchk((/1,ni/),(/1,nj/),angle,angchk)
  angchk(1,:) = -angchk(ni,:)
  ! reverse angle for MOM6
  angchk = -angchk
  write(logmsg,'(a,2f12.2)')'ANGCHK min,max: ',minval(angchk),maxval(angchk)
  print '(a)',trim(logmsg)
  write(logmsg,'(a,2f12.2)')'ANGCHK edges i=1,i=ni: ',angchk(1,nj),angchk(ni,nj)
  print '(a)',trim(logmsg)
 
  !---------------------------------------------------------------------
  ! For the 1/4deg grid, hte at j=720 and j = 1440 is identically=0.0 for
  ! j > 840 (64.0N). These are land points, but since CICE uses hte to
  ! generate remaining variables, setting them to zero will cause problems
  ! For 1deg grid, hte at ni/2 and ni are very small O~10-12, so test for
  ! hte < 1.0
  !---------------------------------------------------------------------
 
  write(logmsg,'(a,2e12.5)')'min vals of hte at folds ', minval(hte(ni/2,:)),minval(hte(ni,:))
  print '(a)',trim(logmsg)
  do j = 1,nj
     ii = ni/2
     if(hte(ii,j) .le. 1.0)hte(ii,j) = 0.5*(hte(ii-1,j) + hte(ii+1,j))
     ii = ni
     if(hte(ii,j) .le. 1.0)hte(ii,j) = 0.5*(hte(ii-1,j) + hte(   1,j))
  enddo
  write(logmsg,'(a,2e12.5)')'min vals of hte at folds ', minval(hte(ni/2,:)),minval(hte(ni,:))
  print '(a)',trim(logmsg)
 
  !---------------------------------------------------------------------
  ! find required extended values for setting all vertices
  !---------------------------------------------------------------------
 
  if(debug)call checkseam
 
  do i = 1,ni
     i2 = ipole(2)+(ipole(1)-i)+1
     xlonct(i) = lonct(i2,nj)
     xlatct(i) = latct(i2,nj)
  enddo
 
  do i = 1,ni
     i2 = ipole(2)+(ipole(1)-i)
     if(i2 .lt. 1)i2 = ni
     xloncu(i) = loncu(i2,nj)
     xlatcu(i) = latcu(i2,nj)
  enddo
 
  if(debug)call checkxlatlon
 
  !approx lat at grid bottom
  do i = 1,ni
     dlatbu(i) = latbu(i,1) + 2.0*(latcu(i,1) - latbu(i,1))
     dlatcv(i) = latct(i,1) + 2.0*(latct(i,1) - latcv(i,1))
  enddo
 
  !---------------------------------------------------------------------
  ! fill grid vertices variables
  !---------------------------------------------------------------------
 
  !Ct and Cu grids align in j
  call fill_vertices(2,nj  , ivertct,jvertct, latbu,lonbu, latct_vert,lonct_vert)
  call           fill_bottom(ivertct,jvertct, latbu,lonbu, latct_vert,lonct_vert,dlatbu)
 
  call fill_vertices(2,nj  , ivertcu,jvertcu, latcv,loncv, latcu_vert,loncu_vert)
  call           fill_bottom(ivertcu,jvertcu, latcv,loncv, latcu_vert,loncu_vert,dlatcv)
 
  !Cv and Bu grids align in j
  call fill_vertices(1,nj-1, ivertcv,jvertcv, latcu,loncu, latcv_vert,loncv_vert)
  call              fill_top(ivertcv,jvertcv, latcu,loncu, latcv_vert,loncv_vert, xlatcu, xloncu)
 
  call fill_vertices(1,nj-1, ivertbu,jvertbu, latct,lonct, latbu_vert,lonbu_vert)
  call              fill_top(ivertbu,jvertbu, latct,lonct, latbu_vert,lonbu_vert, xlatct, xlonct)
 
  if(debug)call checkpoint
 
  if(minval(latct_vert) .lt. -1.e3)stop
  if(minval(lonct_vert) .lt. -1.e3)stop
  if(minval(latcu_vert) .lt. -1.e3)stop
  if(minval(loncu_vert) .lt. -1.e3)stop
  if(minval(latcv_vert) .lt. -1.e3)stop
  if(minval(loncv_vert) .lt. -1.e3)stop
  if(minval(latbu_vert) .lt. -1.e3)stop
  if(minval(lonbu_vert) .lt. -1.e3)stop
  deallocate(xlonct, xlatct, xloncu, xlatcu, dlatbu, dlatcv)
 
  !---------------------------------------------------------------------
  ! write out grid file files
  !---------------------------------------------------------------------
 
  ! create a history attribute
  call date_and_time(date=cdate)
  history = 'created on '//trim(cdate)//' from '//trim(fsrc)
 
  ! write fix grid
  fdst = trim(dirout)//'/'//'tripole.mx'//trim(res)//'.nc'
  call write_tripolegrid(trim(fdst))
 
  ! write cice grid
  fdst = trim(dirout)//'/'//'grid_cice_NEMS_mx'//trim(res)//'.nc'
  call write_cicegrid(trim(fdst))
  deallocate(ulon, ulat, htn, hte)
 
  ! write SCRIP files for generation of positional weights
  do k = 1,nv
     cstagger = trim(staggerlocs(k))
     fdst = trim(dirout)//'/'//trim(cstagger)//'.mx'//trim(res)//'_SCRIP.nc'
     logmsg = 'creating SCRIP file '//trim(fdst)
     print '(a)',trim(logmsg)
     call write_scripgrid(trim(fdst),trim(cstagger))
  end do
  deallocate(latcv_vert, loncv_vert)
  deallocate(latcu_vert, loncu_vert)
  deallocate(latbu_vert, lonbu_vert)
 
  ! write SCRIP file with land mask, used for mapped ocean mask
  ! and  mesh creation
  cstagger = trim(staggerlocs(1))
  fdst= trim(dirout)//'/'//trim(cstagger)//'.mx'//trim(res)//'_SCRIP_land.nc'
  logmsg = 'creating SCRIP file '//trim(fdst)
  print '(a)',trim(logmsg)
  call write_scripgrid(trim(fdst),trim(cstagger),imask=int(wet4))
  deallocate(latct_vert, lonct_vert)
 
  !---------------------------------------------------------------------
  ! write lat,lon,depth and mask arrays required by ww3 in creating
  ! mod_def file
  !---------------------------------------------------------------------
 
  write(cnx,i4fmt)nx
  write(form1,'(a)')'('//trim(cnx)//'f14.8)'
  write(form2,'(a)')'('//trim(cnx)//'i2)'
 
  allocate(ww3mask(1:ni,1:nj)); ww3mask = int(wet4)
  allocate(ww3dpth(1:ni,1:nj)); ww3dpth = dp4
 
  where(latct .ge. maximum_lat)ww3mask = 3
  !close last row
  ww3mask(:,nj) = 3
 
  open(unit=21,file=trim(dirout)//'/'//'ww3.mx'//trim(res)//'_x.inp',form='formatted')
  open(unit=22,file=trim(dirout)//'/'//'ww3.mx'//trim(res)//'_y.inp',form='formatted')
  open(unit=23,file=trim(dirout)//'/'//'ww3.mx'//trim(res)//'_bottom.inp',form='formatted')
  open(unit=24,file=trim(dirout)//'/'//'ww3.mx'//trim(res)//'_mapsta.inp',form='formatted')
  ! cice0 .ne. cicen requires obstruction map, should be initialized as zeros (w3grid,ln3032)
  open(unit=25,file=trim(dirout)//'/'//'ww3.mx'//trim(res)//'_obstr.inp',form='formatted')
 
  do j = 1,nj
     write( 21,trim(form1))lonct(:,j)
     write( 22,trim(form1))latct(:,j)
  end do
  do j = 1,nj
     write( 23,trim(form1))ww3dpth(:,j)
     write( 24,trim(form2))ww3mask(:,j)
     !'obsx' and 'obsy' arrays ???
     write( 25,trim(form2))ww3mask(:,j)*0
     write( 25,trim(form2))ww3mask(:,j)*0
  end do
 
  close(21); close(22); close(23); close(24); close(25)
  deallocate(ww3mask); deallocate(ww3dpth)
  deallocate(wet4, wet8)
 
  !---------------------------------------------------------------------
  ! use ESMF regridding to produce mapped ocean mask; first generate
  ! conservative regrid weights from ocean to tiles; then generate the
  ! tiled files containing the mapped ocean mask
  !---------------------------------------------------------------------
 
  do n = 1,size(catm)
     npx = catm(n)
     if (npx < 100) then
        write(atmres,'(a,i2)')'C',npx
     elseif (npx < 1000) then
        write(atmres,'(a,i3)')'C',npx
     else
        write(atmres,'(a,i4)')'C',npx
     end if
 
     method=esmf_regridmethod_conserve
     fsrc = trim(dirout)//'/'//'Ct.mx'//trim(res)//'_SCRIP_land.nc'
     fdst = trim(fv3dir)//'/'//trim(atmres)//'/'//trim(atmres)//'_mosaic.nc'
     fwgt = trim(dirout)//'/'//'Ct.mx'//trim(res)//'.to.'//trim(atmres)//'.nc'
     logmsg = 'creating weight file '//trim(fwgt)
     print '(a)',trim(logmsg)
 
     call esmf_regridweightgen(srcfile=trim(fsrc),dstfile=trim(fdst),         &
          weightfile=trim(fwgt), regridmethod=method,                         &
          unmappedaction=esmf_unmappedaction_ignore, ignoredegenerate=.true., &
          netcdf4fileflag=.true., tilefilepath=trim(fv3dir)//'/'//trim(atmres)//'/', rc=rc)
     if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
          line=__line__, file=__file__)) call esmf_finalize(endflag=esmf_end_abort)
 
     logmsg = 'creating mapped ocean mask for '//trim(atmres)
     print '(a)',trim(logmsg)
     call make_frac_land(trim(fsrc), trim(fwgt))
  end do
 
  !---------------------------------------------------------------------
  ! use ESMF to create positional weights for mapping a field from its
  ! native stagger location (Cu,Cv,Bu) onto the center (Ct) grid location
  !---------------------------------------------------------------------
 
  method=esmf_regridmethod_bilinear
  fdst = trim(dirout)//'/'//'Ct.mx'//trim(res)//'_SCRIP.nc'
  do k = 2,nv
     cstagger = trim(staggerlocs(k))
     fsrc = trim(dirout)//'/'//trim(cstagger)//'.mx'//trim(res)//'_SCRIP.nc'
     fwgt = trim(dirout)//'/'//'tripole.mx'//trim(res)//'.'//trim(cstagger)//'.to.Ct.bilinear.nc'
     logmsg = 'creating weight file '//trim(fwgt)
     print '(a)',trim(logmsg)
 
     call esmf_regridweightgen(srcfile=trim(fsrc),dstfile=trim(fdst), &
          weightfile=trim(fwgt), regridmethod=method,                 &
          ignoredegenerate=.true.,                                    &
          unmappedaction=esmf_unmappedaction_ignore, rc=rc)
     if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
          line=__line__, file=__file__)) call esmf_finalize(endflag=esmf_end_abort)
  end do
 
  !---------------------------------------------------------------------
  ! use ESMF to create positional weights for mapping a field from the
  ! center (Ct) grid location back to the native stagger location
  ! (Cu,Cv,Bu).
  !---------------------------------------------------------------------
 
  method=esmf_regridmethod_bilinear
  fsrc = trim(dirout)//'/'//'Ct.mx'//trim(res)//'_SCRIP.nc'
  do k = 2,nv
     cstagger = trim(staggerlocs(k))
     fdst = trim(dirout)//'/'//trim(cstagger)//'.mx'//trim(res)//'_SCRIP.nc'
     fwgt = trim(dirout)//'/'//'tripole.mx'//trim(res)//'.Ct.to.'//trim(cstagger)//'.bilinear.nc'
     logmsg = 'creating weight file '//trim(fwgt)
     print '(a)',trim(logmsg)
 
     call esmf_regridweightgen(srcfile=trim(fsrc),dstfile=trim(fdst), &
          weightfile=trim(fwgt), regridmethod=method,                 &
          ignoredegenerate=.true.,                                    &
          unmappedaction=esmf_unmappedaction_ignore, rc=rc)
     if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
          line=__line__, file=__file__)) call esmf_finalize(endflag=esmf_end_abort)
  end do
 
  !---------------------------------------------------------------------
  !
  !---------------------------------------------------------------------
 
  if(do_postwgts)call make_postwgts
 
  !---------------------------------------------------------------------
  ! clean up
  !---------------------------------------------------------------------
 
  deallocate(x, y, dx, dy)
  deallocate(areact, anglet, angle, angchk)
  deallocate(latct, lonct)
  deallocate(latcv, loncv)
  deallocate(latcu, loncu)
  deallocate(latbu, lonbu)
 
end program gen_fixgrid