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_angq, find_ang
  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, dirsrc, dirout, atmres, 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,i2,j2
  integer :: ii,jj
  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 resolution '//trim(atmres)
  print '(a,i6)',' fv3 tile grid size ',npx
  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)

  !print *,minval(wet8),maxval(wet8)
  !print *,minval(wet4),maxval(wet4)

  !---------------------------------------------------------------------
  ! 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)
  !print *,'super grid size ',size(y,1),size(y,2)
  !print *,'max lat in super grid ',maxval(y)
  sg_maxlat = maxval(y)

  !---------------------------------------------------------------------
  ! find the angle on corners---this requires the supergrid
  !---------------------------------------------------------------------

  call find_angq

  !---------------------------------------------------------------------
  ! 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
        !in rad already
        angle(i,j) = -angq(i2,j2)
        !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

  !---------------------------------------------------------------------
  ! find the angle on centers---this does not requires the supergrid
  !---------------------------------------------------------------------

  call find_ang
  print *,'ANGLET ',minval(anglet),maxval(anglet)
  print *,'ANGLE  ',minval(angle),maxval(angle)

  !---------------------------------------------------------------------
  ! 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)

  !---------------------------------------------------------------------
  !
  !---------------------------------------------------------------------

  where(lonct .lt. 0.0)lonct = lonct + 360._dbl_kind
  where(loncu .lt. 0.0)loncu = loncu + 360._dbl_kind
  where(loncv .lt. 0.0)loncv = loncv + 360._dbl_kind
  where(lonbu .lt. 0.0)lonbu = lonbu + 360._dbl_kind

  !---------------------------------------------------------------------
  ! some basic error checking
  ! find the i-th index of the poles at 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)

  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 grids; only the Ct is required, the remaining
  ! staggers are used only in the postweights generation
  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 = 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
  !---------------------------------------------------------------------

  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))

  !---------------------------------------------------------------------
  ! use ESMF to find the tripole:tripole weights for creation
  ! of CICE ICs; the source grid is always mx025; don't create this
  ! file if destination is also mx025
  !---------------------------------------------------------------------

  if(trim(res) .ne. '025') then
     fsrc = trim(dirout)//'/'//'Ct.mx025_SCRIP.nc'
     inquire(file=trim(fsrc), exist=fexist)
     if (fexist ) then
        method=esmf_regridmethod_nearest_stod
        fdst = trim(dirout)//'/'//'Ct.mx'//trim(res)//'_SCRIP.nc'
        fwgt = trim(dirout)//'/'//'tripole.mx025.Ct.to.mx'//trim(res)//'.Ct.neareststod.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)
     else
        logmsg = 'ERROR: '//trim(fsrc)//' is required to generate tripole:triple weights'
        print '(a)',trim(logmsg)
        stop
     end if
  end if

  !---------------------------------------------------------------------
  !
  !---------------------------------------------------------------------

  if(do_postwgts)call make_postwgts

  !---------------------------------------------------------------------
  ! clean up
  !---------------------------------------------------------------------

  deallocate(x,y, angq, dx, dy, xsgp1, ysgp1)
  deallocate(areact, anglet, angle)
  deallocate(latct, lonct)
  deallocate(latcv, loncv)
  deallocate(latcu, loncu)
  deallocate(latbu, lonbu)

end program gen_fixgrid