add make_topo scripts to generate BCs topography files

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt

@@ -38,6 +38,8 @@ else ()
   target_compile_definitions(${this} PRIVATE STDC)
 endif ()
 
+esma_add_subdirectories (@ncar_topo utils_topo)
+
 ecbuild_add_executable (TARGET CombineRasters.x SOURCES CombineRasters.F90 LIBS MAPL ${this})
 ecbuild_add_executable (TARGET mkCatchParam.x SOURCES mkCatchParam.F90 LIBS MAPL ${this} OpenMP::OpenMP_Fortran)
 ecbuild_add_executable (TARGET mkCubeFVRaster.x SOURCES mkCubeFVRaster.F90 LIBS MAPL ${this})

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/CMakeLists.txt

@@ -0,0 +1,31 @@
+#--------------------
+# Copy include files that are used by other libraries.
+# We could leave these in the source directory, and just broaden the search path
+# in the other libaries, but this make it explicit which aspects are externally
+# used.
+
+ecbuild_add_executable (TARGET generate_scrip_cube_topo.x SOURCES generate_scrip_cube.F90 geompack.F90)
+target_link_libraries (generate_scrip_cube_topo.x PRIVATE MPI::MPI_Fortran esmf)
+# CMake has an OpenMP issue with NAG Fortran: https://gitlab.kitware.com/cmake/cmake/-/issues/21280
+if (NOT CMAKE_Fortran_COMPILER_ID MATCHES "NAG")
+  target_link_libraries(generate_scrip_cube_topo.x PRIVATE OpenMP::OpenMP_Fortran)
+endif ()
+
+ecbuild_add_executable (TARGET convert_bin_to_netcdf_topo.x SOURCES convert_bin_to_netcdf.F90)
+target_link_libraries (convert_bin_to_netcdf_topo.x PRIVATE MPI::MPI_Fortran NetCDF::NetCDF_Fortran)
+# CMake has an OpenMP issue with NAG Fortran: https://gitlab.kitware.com/cmake/cmake/-/issues/21280
+if (NOT CMAKE_Fortran_COMPILER_ID MATCHES "NAG")
+  target_link_libraries(convert_bin_to_netcdf_topo.x PRIVATE OpenMP::OpenMP_Fortran)
+endif ()
+
+ecbuild_add_executable (TARGET convert_to_gmao_output_topo.x SOURCES convert_to_gmao_output.F90)
+target_link_libraries (convert_to_gmao_output_topo.x PRIVATE MPI::MPI_Fortran NetCDF::NetCDF_Fortran)
+# CMake has an OpenMP issue with NAG Fortran: https://gitlab.kitware.com/cmake/cmake/-/issues/21280
+if (NOT CMAKE_Fortran_COMPILER_ID MATCHES "NAG")
+  target_link_libraries(convert_to_gmao_output_topo.x PRIVATE OpenMP::OpenMP_Fortran)
+endif ()
+
+install(PROGRAMS scrip_to_cube_topo.py DESTINATION bin)
+install(PROGRAMS scrip_to_restart_topo.py DESTINATION bin)
+install(PROGRAMS generate_topo.sh DESTINATION bin)
+install(PROGRAMS make_topo.py DESTINATION bin)

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/convert_bin_to_netcdf.F90

@@ -0,0 +1,128 @@
+program create_example
+use netcdf
+use, intrinsic :: iso_fortran_env, only: REAL64
+implicit none
+
+
+character(len=512) :: fin,fout,str,fncar
+integer :: im_world,jm_world
+integer :: varid, lonid, latid
+integer :: i, j, nc,xid,yid
+integer :: zid,gwdid,trbid
+integer :: ncid,rc
+integer :: dimids(2)
+integer :: status
+
+integer :: nargs
+
+logical :: doNcar,doGEOS
+
+integer :: ntiles
+real, allocatable :: z1d(:)
+real(REAL64), allocatable :: xdim(:),ydim(:)
+real, allocatable :: a(:,:)
+logical :: isCube
+
+nargs = command_argument_count()
+
+doNCAR=.false.
+doGEOS=.false.
+isCube = .true.
+do i=1,nargs
+     call get_command_argument(i,str)
+     select case(trim(str))
+     case ('-i','--input')
+         call get_command_argument(i+1,fin)
+     case ('-o','--output')
+         call get_command_argument(i+1,fout)
+         doGEOS=.true.
+     case ('--im')
+          call get_command_argument(i+1,str)
+          read(str,'(I10)')im_world
+     case ('--jm')
+          call get_command_argument(i+1,str)
+          read(str,'(I10)')jm_world
+          isCube = .false.
+     case ('--ncar')
+          call get_command_argument(i+1,fncar)
+          doNCAR=.true.
+     end select
+enddo
+
+if (isCube) jm_world = im_world*6
+
+allocate(a(im_world,jm_world))
+open(file=fin,unit=21,form='unformatted')
+read(21)a
+close(21)
+
+if (doGEOS) then
+
+   call check( nf90_create(fout, NF90_NETCDF4,ncid),"error")
+   call check( nf90_def_dim(ncid,"Xdim",im_world,lonid),"error")
+   call check( nf90_def_var(ncid,"Xdim",NF90_DOUBLE,(/lonid/),xid),"error")
+   call check( nf90_put_att(ncid,xid,"units","degrees_east"),"error")
+   call check( nf90_def_dim(ncid,"Ydim",jm_world,latid),"error")
+   call check( nf90_def_var(ncid,"Ydim",NF90_DOUBLE,(/latid/),yid),"error")
+   call check( nf90_put_att(ncid,yid,"units","degrees_north"),"error")
+   call check( nf90_def_var(ncid,"z",NF90_FLOAT,(/lonid,latid/),varid),"error")
+   call check( nf90_put_att(ncid,varid,"units","m"),"error")
+   call check( nf90_put_att(ncid,varid,"long_name","height above sea level"),"error")
+
+   call check( nf90_enddef(ncid),"error")
+
+   allocate(xdim(im_world),ydim(jm_world))
+   do i=1,im_world
+   xdim(i)=i
+   enddo
+   do j=1,jm_world
+   ydim(j)=j
+   enddo
+
+   call check(nf90_put_var(ncid,xid,xdim),"error")
+   call check(nf90_put_var(ncid,yid,ydim),"error")
+   call check(nf90_put_var(ncid,varid,a),"error")
+   call check(nf90_close(ncid),"error")
+
+end if
+
+if (doNCAR) then
+
+   ntiles=im_world*jm_world
+   allocate(z1d(ntiles))
+   call check( nf90_create(fncar, NF90_NETCDF4,ncid),"error")
+   call check( nf90_def_dim(ncid,"ncol",ntiles,xid),"error")
+   call check( nf90_def_var(ncid,"PHIS",NF90_DOUBLE,(/xid/),varid),"error")
+   call check( nf90_put_att(ncid,varid,"long_name","height"),"error")
+   call check( nf90_put_att(ncid,varid,"units","m"),"error")
+   call check( nf90_enddef(ncid),"error")
+
+   nc=0
+   do j=1,jm_world
+   do i=1,im_world
+      nc=nc+1
+      z1d(nc)=a(i,j)
+   enddo
+   enddo
+
+   call check(nf90_put_var(ncid,varid,z1d),"error")
+
+end if
+
+contains
+
+subroutine check(status,loc)
+
+      integer, intent ( in) :: status
+      character(len=*), intent ( in) :: loc
+
+      if(status /= NF90_noerr) then
+         write (*,*) "Error at ", loc
+         write (*,*) nf90_strerror(status)
+         stop "Stopped"
+      end if
+
+end subroutine check
+
+end program create_example
+

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/utils_topo/convert_to_gmao_output.F90

@@ -0,0 +1,209 @@
+program create_example
+use netcdf
+implicit none
+
+
+character(len=512) :: fin,str,fonc4,fobin,fsnc4,fsbin
+integer :: im_world,jm_world
+integer :: varid, lonid, latid
+integer :: i, j, nc
+integer :: zid,gwdid,trbid
+integer :: ncid,rc
+integer :: dimids(2)
+integer :: status
+
+real :: g
+real, allocatable :: z1d(:),turb1d(:),gwd1d(:)
+
+real, allocatable :: a(:,:)
+real  :: asum
+
+logical :: isLL
+integer :: nargs
+
+nargs = command_argument_count()
+
+isLL = .false.
+do i=1,nargs
+     call get_command_argument(i,str)
+     select case(trim(str))
+     case ('-i','--input')
+         call get_command_argument(i+1,fin)
+     case ('--im')
+          call get_command_argument(i+1,str)
+          read(str,'(I10)')im_world
+          jm_world=im_world*6
+     case ('--jm')
+          call get_command_argument(i+1,str)
+          read(str,'(I10)')jm_world
+          isLL=.true.
+     end select
+enddo
+
+g = 9.80616
+
+allocate(z1d(im_world*jm_world),turb1d(im_world*jm_world),gwd1d(im_world*jm_world))
+allocate(a(im_world,jm_world))
+
+call check(nf90_open(fin,NF90_NOWRITE,ncid),"open")
+call check(nf90_inq_varid(ncid,'PHIS',varid),"find phis")
+call check(nf90_get_var(ncid,varid,z1d),"read phis")
+call check(nf90_inq_varid(ncid,'SGH',varid),"find sgh")
+call check(nf90_get_var(ncid,varid,gwd1d),"read sgh")
+call check(nf90_inq_varid(ncid,'SGH30',varid),"find sgh30")
+call check(nf90_get_var(ncid,varid,turb1d),"read sgh30")
+call check(nf90_close(ncid),"close")
+
+if (isLL) then
+   if (im_world < 288) then
+      write(fsnc4,"(i3.3,'x',i2.2,'.nc4')")im_world,jm_world
+      write(fsbin,"(i3.3,'x',i2.2,'.data')")im_world,jm_world
+   else if (im_world < 1152) then
+      write(fsnc4,"(i3.3,'x',i3.3,'.nc4')")im_world,jm_world
+      write(fsbin,"(i3.3,'x',i3.3,'.data')")im_world,jm_world
+   else
+      write(fsnc4,"(i4.4,'x',i3.3,'.nc4')")im_world,jm_world
+      write(fsbin,"(i4.4,'x',i3.3,'.data')")im_world,jm_world
+   end if
+else
+   if (im_world < 168) then
+      write(fsnc4,"(i2.2,'x',i3.3,'.nc4')")im_world,jm_world
+      write(fsbin,"(i2.2,'x',i3.3,'.data')")im_world,jm_world
+   else if (im_world >= 168 .and. im_world < 1000) then
+      write(fsnc4,"(i3.3,'x',i4.4,'.nc4')")im_world,jm_world
+      write(fsbin,"(i3.3,'x',i4.4,'.data')")im_world,jm_world
+   else if (im_world >= 1000 .and. im_world < 1666) then
+      write(fsnc4,"(i4.4,'x',i4.4,'.nc4')")im_world,jm_world
+      write(fsbin,"(i4.4,'x',i4.4,'.data')")im_world,jm_world
+   else if (im_world >= 1666) then
+      write(fsnc4,"(i4.4,'x',i5.5,'.nc4')")im_world,jm_world
+      write(fsbin,"(i4.4,'x',i5.5,'.data')")im_world,jm_world
+   endif
+end if
+
+
+! mean height
+fonc4="gmted_DYN_ave_"//trim(fsnc4)
+fobin="gmted_DYN_ave_"//trim(fsbin)
+call check(nf90_create(fonc4, IOR(NF90_CLOBBER,NF90_NETCDF4),ncid),"error")
+open(file=fobin,unit=21,form='unformatted')
+call check(nf90_def_dim(ncid,"lon",im_world,lonid),"error")
+call check(nf90_def_dim(ncid,"lat",jm_world,latid),"error")
+call check(nf90_def_var(ncid,"z",NF90_FLOAT,(/lonid,latid/),varid),"error")
+call check(nf90_enddef(ncid),"error")
+
+nc=0
+do j=1,jm_world
+do i=1,im_world
+   nc=nc+1
+   !ncar program outputs PHIS, divide by g
+   a(i,j)=z1d(nc)/g
+enddo
+enddo
+
+if (isLL) then
+   asum=0.0
+   do i=1,im_world
+      asum=asum+a(i,1)
+   enddo
+   a(:,1)=asum/float(im_world)
+   asum=0.0
+   do i=1,im_world
+      asum=asum+a(i,jm_world)
+   enddo
+   a(:,jm_world)=asum/float(im_world)
+end if
+
+call check(nf90_put_var(ncid,varid,a),"error")
+call check(nf90_close(ncid),"error")
+
+write(21)a
+close(21)
+
+! GWD
+fonc4="gmted_GWD_var_"//trim(fsnc4)
+fobin="gmted_GWD_var_"//trim(fsbin)
+call check(nf90_create(fonc4, IOR(NF90_CLOBBER,NF90_NETCDF4),ncid),"error")
+open(file=fobin,unit=21,form='unformatted')
+call check(nf90_def_dim(ncid,"lon",im_world,lonid),"error")
+call check(nf90_def_dim(ncid,"lat",jm_world,latid),"error")
+call check(nf90_def_var(ncid,"gwd",NF90_FLOAT,(/lonid,latid/),varid),"error")
+call check(nf90_enddef(ncid),"error")
+nc=0
+do j=1,jm_world
+do i=1,im_world
+   nc=nc+1
+   a(i,j)=gwd1d(nc)
+enddo
+enddo
+if (isLL) then
+   asum=0.0
+   do i=1,im_world
+      asum=asum+a(i,1)
+   enddo
+   a(:,1)=asum/float(im_world)
+   asum=0.0
+   do i=1,im_world
+      asum=asum+a(i,jm_world)
+   enddo
+   a(:,jm_world)=asum/float(im_world)
+end if
+
+call check(nf90_put_var(ncid,varid,a),"error")
+call check(nf90_close(ncid),"error")
+
+write(21)a
+close(21)
+
+! TURB
+fonc4="gmted_TRB_var_"//trim(fsnc4)
+fobin="gmted_TRB_var_"//trim(fsbin)
+call check(nf90_create(fonc4, IOR(NF90_CLOBBER,NF90_NETCDF4),ncid),"error")
+open(file=fobin,unit=21,form='unformatted')
+call check(nf90_def_dim(ncid,"lon",im_world,lonid),"error")
+call check(nf90_def_dim(ncid,"lat",jm_world,latid),"error")
+call check(nf90_def_var(ncid,"trb",NF90_FLOAT,(/lonid,latid/),varid),"error")
+call check(nf90_enddef(ncid),"error")
+nc=0
+do j=1,jm_world
+do i=1,im_world
+   nc=nc+1
+   a(i,j)=turb1d(nc)
+enddo
+enddo
+if (isLL) then
+   asum=0.0
+   do i=1,im_world
+      asum=asum+a(i,1)
+   enddo
+   a(:,1)=asum/float(im_world)
+   asum=0.0
+   do i=1,im_world
+      asum=asum+a(i,jm_world)
+   enddo
+   a(:,jm_world)=asum/float(im_world)
+end if
+
+call check(nf90_put_var(ncid,varid,a),"error")
+call check(nf90_close(ncid),"error")
+
+write(21)a
+close(21)
+
+contains
+
+subroutine check(status,loc)
+
+      integer, intent ( in) :: status
+      character(len=*), intent ( in) :: loc
+
+      if(status /= NF90_noerr) then
+         write (*,*) "Error at ", loc
+         write (*,*) nf90_strerror(status)
+         stop "Stopped"
+      end if
+
+end subroutine check
+
+end program create_example
+