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)
 mastertask = .false.
 if (localpet == 0) mastertask=.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')

  if(mastertask) then
    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 *
  end if

  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

  if(mastertask) then
    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 *
  end if

  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)
  if(mastertask) then
    print '(a)', 'reading ocean mask from '//trim(fsrc)
    if(rc .ne. 0)print '(a)', 'nf90_open = '//trim(nf90_strerror(rc))
  end if

  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)
  if(mastertask) then
    print '(a)', 'reading ocean topography from '//trim(fsrc)
    if(rc .ne. 0)print '(a)', 'nf90_open = '//trim(nf90_strerror(rc))
  end if

  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)

  !print *,minval(dp8),maxval(dp8)
  !print *,minval(dp4),maxval(dp4)

  if(editmask)then
!---------------------------------------------------------------------
!  apply topoedits run time mask changes if required for this config
!---------------------------------------------------------------------

   if(trim(editsfile)  == 'none')then
    print '(a)', 'Need a valid editsfile to make mask edits '
    stop
   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
!---------------------------------------------------------------------

   fsrc = trim(dirsrc)//'/'//trim(editsfile)
   if(editmask)fsrc = trim(dirout)//'/'//'ufs.'//trim(editsfile)
   call apply_topoedits(fsrc)

!---------------------------------------------------------------------
! read MOM6 supergrid file
!---------------------------------------------------------------------

  fsrc = trim(dirsrc)//'/'//'ocean_hgrid.nc'

  rc = nf90_open(fsrc, nf90_nowrite, ncid)
  if(mastertask) then
    print '(a)', 'reading supergrid from '//trim(fsrc)
    if(rc .ne. 0)print '(a)', 'nf90_open = '//trim(nf90_strerror(rc))
  end if
  
  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

    if(mastertask) then
      print *,'ANGLET ',minval(anglet),maxval(anglet)
      print *,'ANGLE  ',minval(angle),maxval(angle)
    end if

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

   if(mastertask) then
     write(logmsg,'(a,2e12.5)')'min vals of hte at folds ', &
                              minval(hte(ni/2,:)),minval(hte(ni,:))
     print '(a)',trim(logmsg)
   end if
   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
   if(mastertask) then
     write(logmsg,'(a,2e12.5)')'min vals of hte at folds ', &
                              minval(hte(ni/2,:)),minval(hte(ni,:))
     print '(a)',trim(logmsg)
   end if

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

  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
    if(mastertask) then
      write(logmsg,'(a,2i6,2f12.2)')'poles found at i = ',ipole,latbu(ipole(1),nj), &
                                                                latbu(ipole(2),nj)
      print '(a)',trim(logmsg)
    end if

    if(mastertask .and. 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(mastertask .and. 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(mastertask .and. 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'
    if(mastertask) then
      logmsg = 'creating SCRIP file '//trim(fdst)
      print '(a)',trim(logmsg)
    end if
    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'
   if(mastertask) then
     logmsg = 'creating SCRIP file '//trim(fdst)
     print '(a)',trim(logmsg)
   end if
   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'
   if(mastertask) then
     logmsg = 'creating weight file '//trim(fwgt)
     print '(a)',trim(logmsg)
   end if

   call esmf_regridweightgen(srcfile=trim(fsrc),dstfile=trim(fdst), &
                        weightfile=trim(fwgt), regridmethod=method, &
                        unmappedaction=esmf_unmappedaction_ignore, &
                        ignoredegenerate=.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)

   if(mastertask) then
     logmsg = 'creating mapped ocean mask for '//trim(atmres)
     print '(a)',trim(logmsg)
   end if
   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'
       if(mastertask) then
         logmsg = 'creating weight file '//trim(fwgt)
         print '(a)',trim(logmsg)
       end if
       
       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
       if(mastertask) then
         logmsg = 'ERROR: '//trim(fsrc)//' is required to generate tripole:triple weights'
         print '(a)',trim(logmsg)
       end if
       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