First commit to bring coupled gridgen into ufs utils (ufs-community#143)

build_all.sh

@@ -26,8 +26,8 @@ else
  set -x
 fi
 
-CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=../ -DEMC_EXEC_DIR=ON"
-#CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=../ -DEMC_EXEC_DIR=ON -DENABLE_DOCS=ON"
+#CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=../ -DEMC_EXEC_DIR=ON"
+CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=../ -DEMC_EXEC_DIR=ON -DENABLE_DOCS=ON"
 
 rm -fr ./build
 mkdir ./build && cd ./build

modulefiles/build.hera.intel

@@ -28,4 +28,5 @@ module load png/1.6.35
 module load hdf5/1.10.6
 module load netcdf/4.7.4
 module load nccmp/1.8.7.0
-module load esmf/8_1_0_beta_snapshot_27
+#module load esmf/8_1_0_beta_snapshot_27
+module load esmf/8_2_0

sorc/CMakeLists.txt

@@ -18,3 +18,4 @@ add_subdirectory(sfc_climo_gen.fd)
 add_subdirectory(vcoord_gen.fd)
 add_subdirectory(fvcom_tools.fd)
 add_subdirectory(gblevents.fd)
+add_subdirectory(cpld_gridgen.fd)

sorc/cpld_gridgen.fd/CMakeLists.txt

@@ -0,0 +1,61 @@
+# This is the CMake build file for the chgres_cube utility in the
+# UFS_UTILS package.
+#
+# George Gayno, Mark Potts, Kyle Gerheiser
+
+set(lib_src
+    angles.F90
+    charstrings.F90
+    cicegrid.F90
+    debugprint.F90
+    gengrid_kinds.F90
+    grdvars.F90
+    inputnml.F90
+    mapped_mask.F90
+    postwgts.F90
+    scripgrid.F90
+    topoedits.F90
+    tripolegrid.F90
+    vartypedefs.F90
+    vertices.F90)
+
+set(exe_src gen_fixgrid.F90)
+
+if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$")
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -assume byterecl")
+elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$")
+  set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8")
+
+  # Turn on this argument mismatch flag for gfortran10.
+  if(CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 10)
+     set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fallow-argument-mismatch")
+  endif()
+endif()
+
+set(exe_name cpld_gridgen)
+
+add_library(cpld_gridgen_lib STATIC ${lib_src})
+add_executable(${exe_name} ${exe_src})
+
+set(mod_dir "${CMAKE_CURRENT_BINARY_DIR}/mod")
+set_target_properties(cpld_gridgen_lib PROPERTIES Fortran_MODULE_DIRECTORY ${mod_dir})
+target_include_directories(cpld_gridgen_lib INTERFACE ${mod_dir})
+
+target_link_libraries(
+  cpld_gridgen_lib
+  PUBLIC
+  esmf
+  MPI::MPI_Fortran)
+
+if(OpenMP_Fortran_FOUND)
+  target_link_libraries(${exe_name} PUBLIC OpenMP::OpenMP_Fortran)
+endif()
+
+target_link_libraries(${exe_name} PRIVATE cpld_gridgen_lib)
+
+install(TARGETS ${exe_name} RUNTIME DESTINATION ${exec_dir})
+
+# If doxygen documentation we enabled, build it.
+if(ENABLE_DOCS)
+  add_subdirectory(docs)  
+endif()

sorc/cpld_gridgen.fd/angles.F90

@@ -0,0 +1,187 @@
+module angles
+
+  use gengrid_kinds, only : dbl_kind, int_kind
+  use grdvars,       only : ni,nj,nx,ny
+  use grdvars,       only : x,y,xsgp1,ysgp1,sg_maxlat
+  use grdvars,       only : latBu,lonBu,lonCt
+  use grdvars,       only : angq,anglet
+  use grdvars,       only : mastertask, debug
+
+  implicit none
+
+  contains
+
+  subroutine find_angq
+
+    ! local variables
+    integer :: i,j,i1,i2,m,n
+
+    ! pole locations on SG
+    integer(int_kind) :: ipolesg(2)
+
+    ! from geolonB fix in MOM6
+    real(dbl_kind) :: len_lon ! The periodic range of longitudes, usually 360 degrees.
+    real(dbl_kind) :: pi_720deg ! One quarter the conversion factor from degrees to radians.
+    real(dbl_kind) :: lonB(2,2)  ! The longitude of a point, shifted to have about the same value.
+    real(dbl_kind) :: lon_scale = 0.0
+
+    real(dbl_kind) :: modulo_around_point
+!---------------------------------------------------------------------
+! to find angleq on seam, replicate supergrid values across seam
+!---------------------------------------------------------------------
+
+     angq = 0.0
+    xsgp1 = 0.0; ysgp1 = 0.0
+    !pole on supergrid
+    ipolesg = -1
+        j = ny
+    do i = 1,nx/2
+     if(y(i,j) .eq. sg_maxlat)ipolesg(1) = i
+    enddo
+    do i = nx/2+1,nx
+     if(y(i,j) .eq. sg_maxlat)ipolesg(2) = i
+    enddo
+    if(mastertask .and. debug)print *,'poles found at ',ipolesg
+
+    xsgp1(:,0:ny) = x(:,0:ny)
+    ysgp1(:,0:ny) = y(:,0:ny)
+
+    !check
+    do i = ipolesg(1)-5,ipolesg(1)+5
+     i2 = ipolesg(2)+(ipolesg(1)-i)+1
+     if(mastertask .and. debug)print *,i,i2
+    enddo
+     print *
+    do i = ipolesg(2)-5,ipolesg(2)+5
+     i2 = ipolesg(2)+(ipolesg(1)-i)+1
+     if(mastertask .and. debug)print *,i,i2
+    enddo
+
+    !replicate supergrid across pole
+    do i = 1,nx
+      i2 = ipolesg(2)+(ipolesg(1)-i)
+     xsgp1(i,ny+1) = xsgp1(i2,ny)
+     ysgp1(i,ny+1) = ysgp1(i2,ny)
+    enddo
+
+    !check
+    if(mastertask .and. debug)then
+     j = ny+1
+    i1 = ipolesg(1); i2 = ipolesg(2)-(ipolesg(1)-i1)
+      print *,'replicate X across seam on SG'
+      print *,xsgp1(i1-2,j),xsgp1(i2+2,j)
+      print *,xsgp1(i1-1,j),xsgp1(i2+1,j)
+      print *,xsgp1(i1,  j),xsgp1(i2,  j)
+      print *,xsgp1(i1+1,j),xsgp1(i2-1,j)
+      print *,xsgp1(i1+2,j),xsgp1(i2-2,j)
+
+      print *,'replicate Y across seam on SG'
+      print *,ysgp1(i1-2,j),ysgp1(i2+2,j)
+      print *,ysgp1(i1-1,j),ysgp1(i2+1,j)
+      print *,ysgp1(i1,  j),ysgp1(i2,  j)
+      print *,ysgp1(i1+1,j),ysgp1(i2-1,j)
+      print *,ysgp1(i1+2,j),ysgp1(i2-2,j)
+    end if
+
+!---------------------------------------------------------------------
+! rotation angle on supergrid vertices
+! lonB: x(i-1,j-1) has same relationship to x(i,j) on SG as
+!       geolonT(i,j) has to geolonBu(i,j) on the reduced grid
+!---------------------------------------------------------------------
+
+     ! constants as defined in MOM
+     pi_720deg = atan(1.0) / 180.0
+       len_lon = 360.0
+    do j=1,ny ; do i=1,nx-1
+      do n=1,2 ; do m=1,2
+        lonB(m,n) = modulo_around_point(xsgp1(I+m-2,J+n-2), xsgp1(i-1,j-1), len_lon)
+      enddo ; enddo
+      lon_scale    = cos(pi_720deg*(ysgp1(i-1,j-1) + ysgp1(i+1,j-1) + &
+                                    ysgp1(i-1,j+1) + ysgp1(i+1,j+1)) )
+      angq(i,j)    = atan2(lon_scale*((lonB(1,2) - lonB(2,1)) + (lonB(2,2) - lonB(1,1))), &
+                            ysgp1(i-1,j+1) + ysgp1(i+1,j+1) - &
+                            ysgp1(i-1,j-1) - ysgp1(i+1,j-1) )
+    enddo ; enddo
+
+    !check
+    if(mastertask .and. debug) then
+       j = ny
+      i1 = ipolesg(1); i2 = ipolesg(2)-(ipolesg(1)-i1)
+      print *,'angq along seam on SG'
+      print *,angq(i1-2,j),angq(i2+2,j)
+      print *,angq(i1-1,j),angq(i2+1,j)
+      print *,angq(i1,  j),angq(i2,  j)
+      print *,angq(i1+1,j),angq(i2-1,j)
+      print *,angq(i1+2,j),angq(i2-2,j)
+    end if
+
+  end subroutine find_angq
+
+  subroutine find_ang
+
+    ! local variables
+    integer :: i,j,m,n
+    integer :: ii,jj
+
+    ! from geolonB fix in MOM6
+    real(dbl_kind) :: len_lon ! The periodic range of longitudes, usually 360 degrees.
+    real(dbl_kind) :: pi_720deg ! One quarter the conversion factor from degrees to radians.
+    real(dbl_kind) :: lonB(2,2)  ! The longitude of a point, shifted to have about the same value.
+    real(dbl_kind) :: lon_scale = 0.0
+
+    real(dbl_kind) :: modulo_around_point
+!---------------------------------------------------------------------
+! rotation angle for "use_bugs" = false case from MOM6
+! src/initialization/MOM_shared_initialization.F90 but allow for not
+! having halo values
+! note this does not reproduce sin_rot,cos_rot found in MOM6 output
+! differences are ~O 10-6
+!---------------------------------------------------------------------
+
+     anglet = 0.0
+     pi_720deg = atan(1.0) / 180.0
+       len_lon = 360.0
+      do j=1,nj; do i = 1,ni
+        do n=1,2 ; do m=1,2
+                      jj = J+n-2; ii = I+m-2
+         if(jj .eq. 0)jj = 1
+         if(ii .eq. 0)ii = ni
+          lonB(m,n) = modulo_around_point(LonBu(ii,jj), LonCt(i,j), len_lon)
+        !  lonB(m,n) = modulo_around_point(LonBu(I+m-2,J+n-2), LonCt(i,j), len_lon)
+        enddo ; enddo
+                      jj = j-1; ii = i-1
+         if(jj .eq. 0)jj = 1
+         if(ii .eq. 0)ii = ni
+        lon_scale = cos(pi_720deg*((LatBu(ii,jj) + LatBu(I,J)) + &
+                                   (LatBu(I,jj) + LatBu(ii,J)) ) )
+       anglet(i,j) = atan2(lon_scale*((lonB(1,2) - lonB(2,1)) + (lonB(2,2) - lonB(1,1))), &
+                      (LatBu(ii,J) - LatBu(I,jj)) + &
+                      (LatBu(I,J) - LatBu(ii,jj)) )
+
+        !lon_scale = cos(pi_720deg*((LatBu(I-1,J-1) + LatBu(I,J)) + &
+        !                           (LatBu(I,J-1) + LatBu(I-1,J)) ) )
+       !anglet(i,j) = atan2(lon_scale*((lonB(1,2) - lonB(2,1)) + (lonB(2,2) - lonB(1,1))), &
+       !               (LatBu(I-1,J) - LatBu(I,J-1)) + &
+       !               (LatBu(I,J) - LatBu(I-1,J-1)) )
+      enddo ; enddo
+
+  end subroutine find_ang
+end module angles
+
+! -----------------------------------------------------------------------------
+!> Return the modulo value of x in an interval [xc-(Lx/2) xc+(Lx/2)]
+!! If Lx<=0, then it returns x without applying modulo arithmetic.
+function modulo_around_point(x, xc, Lx) result(x_mod)
+  use gengrid_kinds, only : dbl_kind
+
+  real(dbl_kind), intent(in) :: x  !< Value to which to apply modulo arithmetic
+  real(dbl_kind), intent(in) :: xc !< Center of modulo range
+  real(dbl_kind), intent(in) :: Lx !< Modulo range width
+  real(dbl_kind) :: x_mod          !< x shifted by an integer multiple of Lx to be close to xc.
+
+  if (Lx > 0.0) then
+    x_mod = modulo(x - (xc - 0.5*Lx), Lx) + (xc - 0.5*Lx)
+  else
+    x_mod = x
+  endif
+end function modulo_around_point

sorc/cpld_gridgen.fd/charstrings.F90

@@ -0,0 +1,23 @@
+module charstrings
+
+  use gengrid_kinds, only : CL,CM,CS
+
+  implicit none
+
+  character(len=CL) :: dirsrc, dirout, fv3dir
+  character(len=CS) :: res, atmres
+  character(len=CL) :: logmsg
+
+  character(len=CL) :: maskfile  = 'ocean_mask.nc'
+  character(len=CS) :: maskname  = 'mask'
+  character(len=CL) :: editsfile
+
+  character(len=CL) :: topofile
+  character(len=CS) :: toponame  = 'depth'
+
+  character(len=CL) :: history
+  character(len=CS) :: cdate
+
+  character(len= 2), dimension(4) :: staggerlocs = (/'Ct','Cu','Cv','Bu'/)
+
+end module charstrings

sorc/cpld_gridgen.fd/cicegrid.F90

@@ -0,0 +1,86 @@
+module cicegrid
+
+  use grdvars,       only: ni,nj,ulat,ulon,htn,hte,angle,wet4,mastertask
+  use charstrings,   only: history, logmsg
+  use vartypedefs,   only: maxvars, cicevars, cicevars_typedefine
+  use gengrid_kinds, only: CM
+  use netcdf
+
+  implicit none
+  private
+
+  public write_cicegrid
+
+  contains
+
+  subroutine write_cicegrid(fname)
+
+  character(len=*), intent(in) :: fname
+
+  ! local variables
+
+  integer :: ii,id,rc, ncid, dim2(2)
+  integer :: idimid,jdimid
+
+  character(len=2)  :: vtype
+  character(len=CM) :: vname
+  character(len=CM) :: vlong
+  character(len=CM) :: vunit
+
+!---------------------------------------------------------------------
+! create the netcdf file
+!---------------------------------------------------------------------
+
+  ! define the output variables and file name
+  call cicevars_typedefine
+
+  rc = nf90_create(fname, nf90_write, ncid)
+  if(mastertask) then
+    logmsg = '==> writing CICE grid to '//trim(fname)
+    print '(a)', trim(logmsg)
+    if(rc .ne. 0)print '(a)', 'nf90_create = '//trim(nf90_strerror(rc))
+  end if
+
+  rc = nf90_def_dim(ncid, 'ni', ni, idimid)
+  rc = nf90_def_dim(ncid, 'nj', nj, jdimid)
+
+  do ii = 1,maxvars
+   if(len_trim(cicevars(ii)%var_name) .gt. 0)then
+     vname = trim(cicevars(ii)%var_name)
+     vlong = trim(cicevars(ii)%long_name)
+     vunit = trim(cicevars(ii)%unit_name)
+     vtype = trim(cicevars(ii)%var_type)
+
+     dim2(:) =  (/idimid, jdimid/)
+     if(vtype .eq. 'r8')rc = nf90_def_var(ncid, vname, nf90_double, dim2, id)
+     if(vtype .eq. 'r4')rc = nf90_def_var(ncid, vname, nf90_float,  dim2, id)
+     if(vtype .eq. 'i4')rc = nf90_def_var(ncid, vname, nf90_int,    dim2, id)
+     rc = nf90_put_att(ncid, id,     'units', vunit)
+     rc = nf90_put_att(ncid, id, 'long_name', vlong)
+   end if
+  enddo
+   rc = nf90_put_att(ncid, nf90_global, 'history', trim(history))
+   rc = nf90_enddef(ncid)
+
+  rc = nf90_inq_varid(ncid,  'ulon',      id)
+  rc = nf90_put_var(ncid,        id,    ulon)
+
+  rc = nf90_inq_varid(ncid,  'ulat',      id)
+  rc = nf90_put_var(ncid,        id,    ulat)
+
+  rc = nf90_inq_varid(ncid,   'htn',      id)
+  rc = nf90_put_var(ncid,        id,     htn)
+
+  rc = nf90_inq_varid(ncid,   'hte',      id)
+  rc = nf90_put_var(ncid,        id,     hte)
+
+  rc = nf90_inq_varid(ncid,  'angle',     id)
+  rc = nf90_put_var(ncid,         id,  angle)
+
+  rc = nf90_inq_varid(ncid,    'kmt',        id)
+  rc = nf90_put_var(ncid,         id, int(wet4))
+
+  rc = nf90_close(ncid)
+
+  end subroutine write_cicegrid
+end module cicegrid