[wip] homme DA

components/homme/CMakeLists.txt

@@ -28,6 +28,11 @@ PROJECT(HOMME C Fortran)
 
 INCLUDE(FortranCInterface)
 
+#if machine file sets WITH_DA, then all targets are built with DA (for unit testing)
+IF (WITH_DA)
+ADD_DEFINITIONS(-DHOMMEDA)
+ENDIF ()
+
 IF (${HOMME_BUILD_EXECS})
   OPTION(HOMME_ENABLE_TESTING "Whether tests should be enabled" TRUE)
   IF (${HOMME_ENABLE_TESTING})
@@ -119,7 +124,7 @@ OPTION(BUILD_HOMME_SWEQX "Shallow water equations FEM" ON)
 OPTION(BUILD_HOMME_PREQX "Primitive equations FEM" ON)
 OPTION(BUILD_HOMME_THETA "SB81 NH model with theta and Exner pressure" ON)
 OPTION(BUILD_HOMME_THETA_KOKKOS "SB81 NH model with theta and Exner pressure, Kokkos target" OFF)
-OPTION(BUILD_HOMME_PREQX_ACC "Primitive equations FEM with OpenACC" ON)
+OPTION(BUILD_HOMME_PREQX_ACC "Primitive equations FEM with OpenACC" OFF)
 OPTION(BUILD_HOMME_PREQX_KOKKOS "Primitive equations FEM with Kokkos" OFF)
 OPTION(BUILD_HOMME_SWIM  "Shallow water equations implicit" OFF)
 OPTION(BUILD_HOMME_PRIM  "Primitive equations implicit" OFF)
@@ -567,3 +572,5 @@ ADD_CUSTOM_TARGET(clean-all
    COMMAND ${CMAKE_BUILD_TOOL} clean
    COMMAND ${CMAKE_COMMAND} -P "${CMAKE_CURRENT_SOURCE_DIR}/cmake/clean-cmake-files.cmake"
 )
+
+

components/homme/cmake/machineFiles/chrysalis.cmake

@@ -28,6 +28,8 @@ EXECUTE_PROCESS(COMMAND pnetcdf-config --prefix
 )
 SET (PnetCDF_PATH "${PNCCONFIG_OUTPUT}" CACHE STRING "")
 
+SET (WITH_DA FALSE CACHE BOOL "")
+
 EXECUTE_PROCESS(COMMAND nf-config --prefix
   RESULT_VARIABLE NFCONFIG_RESULT
   OUTPUT_VARIABLE NFCONFIG_OUTPUT

components/homme/dcmip_tests/dcmip2016_test1_baroclinic_wave/theta-l/namelist-r100-topo.nl

@@ -0,0 +1,54 @@
+!
+! preqx: namelist for dcmip2016 test1: moist baroclinic wave
+!_______________________________________________________________________
+&ctl_nl
+  nthreads          = -1                        ! use OMP_NUM_THREADS
+  partmethod        = 4                         ! mesh parition method: 4 = space filling curve
+  topology          = "cube"                    ! mesh type: cubed sphere
+  test_case         = "dcmip2016_test1"         ! test identifier
+  sub_case          = 2
+  ne                = 30                        ! number of elements per cube face
+  qsize             = 6                         ! num tracer fields
+  ndays             = 6
+  statefreq         = 72                        ! number of steps between screen dumps
+  restartfreq       = -1                        ! don't write restart files if < 0
+  runtype           = 0                         ! 0 => new run
+  tstep             = 300                       ! largest timestep in seconds
+  theta_advect_form = 1
+  integration       = 'explicit'                ! explicit time integration
+  tstep_type        = 5                         ! use 9 if running NH
+  dt_remap_factor   = 2
+  dt_tracer_factor  = 1
+  hypervis_scaling  =  3.0                      ! turn on RRM HV
+  nu                = 3.4e-8                    ! for all resolutions
+  nu_top            = 0                         ! default = 2.5e5
+  limiter_option    = 10
+  hypervis_order    = 2                         ! 2 = hyperviscosity
+  hypervis_subcycle = 1                         ! 1 = no hyperviz subcycling
+  moisture          = 'dry'
+  theta_hydrostatic_mode = .true.
+  dcmip16_prec_type = -1                          ! 0=kessler physics
+  dcmip16_pbl_type  = -1                         ! 0=reed-jablonowski pbl, -1 = none
+/
+&vert_nl
+  vfile_mid         = "../vcoord/camm-30.ascii"
+  vfile_int         = "../vcoord/cami-30.ascii"
+/
+&analysis_nl
+!  output_prefix     = "r100-dry-"
+  output_dir        = "./movies/"               ! destination dir for netcdf file
+  output_timeunits  = 2,                        ! 0=timesteps, 1=days, 2=hours, 3=seconds
+  output_frequency  = 24                        ! every N hours
+!  output_varnames1  ='T','ps','pnh','geo','u','v','w','omega','Th','Q','Q2','Q3','Q4','Q5','rho','precl','zeta'   ! variables to write to file
+  output_varnames1  ='T','ps','pnh','geo','u','Th','Q','Q2','Q3','Q4','Q5','precl','zeta'   ! variables to write to file
+  interp_type       = 1                         ! 0=native grid, 1=bilinear
+  output_type       ='netcdf'                   ! netcdf or pnetcdf
+  num_io_procs      = 16
+  interp_nlon       = 1024
+  interp_nlat       = 513
+  interp_gridtype   = 1
+/
+&prof_inparm
+  profile_outpe_num   = 100
+  profile_single_file	= .true.
+/

components/homme/dcmip_tests/dcmip2016_test1_baroclinic_wave/theta-l/namelist-r25-topo.nl

@@ -0,0 +1,54 @@
+!
+! preqx: namelist for dcmip2016 test1: moist baroclinic wave
+!_______________________________________________________________________
+&ctl_nl
+  nthreads          = -1                        ! use OMP_NUM_THREADS
+  partmethod        = 4                         ! mesh parition method: 4 = space filling curve
+  topology          = "cube"                    ! mesh type: cubed sphere
+  test_case         = "dcmip2016_test1"         ! test identifier
+  sub_case          = 2
+  ne                = 120                        ! number of elements per cube face
+  qsize             = 6                         ! num tracer fields
+  ndays             = 6
+  statefreq         = 72                        ! number of steps between screen dumps
+  restartfreq       = -1                        ! don't write restart files if < 0
+  runtype           = 0                         ! 0 => new run
+  tstep             = 75                        ! largest timestep in seconds
+  theta_advect_form = 1
+  integration       = 'explicit'                ! explicit time integration
+  tstep_type        = 5                         ! use 9 if running NH
+  dt_remap_factor   = 2
+  dt_tracer_factor  = 1
+  hypervis_scaling  =  3.0                      ! turn on RRM HV
+  nu                = 3.4e-8                    ! for all resolutions
+  nu_top            = 0                         ! default = 2.5e5
+  limiter_option    = 10
+  hypervis_order    = 2                         ! 2 = hyperviscosity
+  hypervis_subcycle = 1                         ! 1 = no hyperviz subcycling
+  moisture          = 'dry'
+  theta_hydrostatic_mode = .true.
+  dcmip16_prec_type = -1                          ! 0=kessler physics
+  dcmip16_pbl_type  = -1                         ! 0=reed-jablonowski pbl, -1 = none
+/
+&vert_nl
+  vfile_mid         = "../vcoord/camm-30.ascii"
+  vfile_int         = "../vcoord/cami-30.ascii"
+/
+&analysis_nl
+!  output_prefix     = "r100-dry-"
+  output_dir        = "./movies/"               ! destination dir for netcdf file
+  output_timeunits  = 2,                        ! 0=timesteps, 1=days, 2=hours, 3=seconds
+  output_frequency  = 24                        ! every N hours
+!  output_varnames1  ='T','ps','pnh','geo','u','v','w','omega','Th','Q','Q2','Q3','Q4','Q5','rho','precl','zeta'   ! variables to write to file
+  output_varnames1  ='T','ps','pnh','geo','u','Th','Q','Q2','Q3','Q4','Q5','precl','zeta'   ! variables to write to file
+  interp_type       = 1                         ! 0=native grid, 1=bilinear
+  output_type       ='netcdf'                   ! netcdf or pnetcdf
+  num_io_procs      = 16
+  interp_nlon       = 4096
+  interp_nlat       = 2049
+  interp_gridtype   = 1
+/
+&prof_inparm
+  profile_outpe_num   = 100
+  profile_single_file	= .true.
+/

components/homme/dcmip_tests/dcmip2016_test1_baroclinic_wave/theta-l/namelist-rrm-topo.nl

@@ -0,0 +1,55 @@
+!
+! preqx: namelist for dcmip2016 test1: moist baroclinic wave
+!_______________________________________________________________________
+&ctl_nl
+  nthreads          = -1                        ! use OMP_NUM_THREADS
+  partmethod        = 4                         ! mesh parition method: 4 = space filling curve
+  topology          = "cube"                    ! mesh type: cubed sphere
+  test_case         = "dcmip2016_test1"         ! test identifier
+  sub_case          = 2
+  ne                = 0                        ! number of elements per cube face
+  mesh_file         = '/lcrc/group/e3sm/data/inputdata/atm/cam/inic/homme/northamericax4v1.g'
+  qsize             = 6                         ! num tracer fields
+  ndays             = 6
+  statefreq         = 72                        ! number of steps between screen dumps
+  restartfreq       = -1                        ! don't write restart files if < 0
+  runtype           = 0                         ! 0 => new run
+  tstep             = 75                        ! largest timestep in seconds
+  theta_advect_form = 1
+  integration       = 'explicit'                ! explicit time integration
+  tstep_type        = 5                         ! use 9 if running NH
+  dt_remap_factor   = 2
+  dt_tracer_factor  = 1
+  hypervis_scaling               =  3.0
+  nu                = 3.4e-8
+  nu_top            = 0                         ! default = 2.5e5
+  limiter_option    = 10
+  hypervis_order    = 2                         ! 2 = hyperviscosity
+  hypervis_subcycle = 2                         ! is 1 stable?
+  moisture          = 'dry'
+  theta_hydrostatic_mode = .true.
+  dcmip16_prec_type = -1                          ! 0=kessler physics
+  dcmip16_pbl_type  = -1                         ! 0=reed-jablonowski pbl, -1 = none
+/
+&vert_nl
+  vfile_mid         = "../vcoord/camm-30.ascii"
+  vfile_int         = "../vcoord/cami-30.ascii"
+/
+&analysis_nl
+!  output_prefix     = "r100-dry-"
+  output_dir        = "./movies/"               ! destination dir for netcdf file
+  output_timeunits  = 2,                        ! 0=timesteps, 1=days, 2=hours, 3=seconds
+  output_frequency  = 24                        ! every N hours
+!  output_varnames1  ='T','ps','pnh','geo','u','v','w','omega','Th','Q','Q2','Q3','Q4','Q5','rho','precl','zeta'   ! variables to write to file
+  output_varnames1  ='T','ps','pnh','geo','u','Th','min_dx','zeta'   ! variables to write to file
+  interp_type       = 1                         ! 0=native grid, 1=bilinear
+  output_type       ='netcdf'                   ! netcdf or pnetcdf
+  num_io_procs      = 16
+  interp_nlon       = 4096
+  interp_nlat       = 2049
+  interp_gridtype   = 1
+/
+&prof_inparm
+  profile_outpe_num   = 100
+  profile_single_file	= .true.
+/

components/homme/src/preqx/derivative_mod.F90

@@ -8,6 +8,7 @@ module derivative_mod
        gradient_sphere_wk_testcov, gradient_sphere_wk_testcontra, ugradv_sphere, vorticity_sphere, vorticity_sphere_diag,    &
        curl_sphere,curl_sphere_wk_testcov, vlaplace_sphere_wk, element_boundary_integral, &
        limiter_optim_iter_full, limiter_clip_and_sum, laplace_sphere_wk, divergence_sphere_wk, gradient_sphere,              &
-       divergence_sphere, laplace_z
+       divergence_sphere, laplace_z, &
+       get_deriv
   implicit none
 end module derivative_mod

components/homme/src/preqx_kokkos/derivative_mod.F90

@@ -8,6 +8,6 @@ module derivative_mod
        gradient_sphere_wk_testcov, gradient_sphere_wk_testcontra, ugradv_sphere, vorticity_sphere, vorticity_sphere_diag,    &
        curl_sphere,curl_sphere_wk_testcov, vlaplace_sphere_wk, element_boundary_integral, &
        limiter_optim_iter_full, limiter_clip_and_sum, laplace_sphere_wk, divergence_sphere_wk, gradient_sphere,              &
-       divergence_sphere, laplace_z
+       divergence_sphere, laplace_z, get_deriv
   implicit none
 end module derivative_mod

components/homme/src/share/cube_mod.F90

@@ -788,11 +788,15 @@ subroutine coriolis_init_atomic(elem)
              if ( rotate_grid /= 0) then
                 lat = elem%spherep(i,j)%lat
                 lon = elem%spherep(i,j)%lon
-             	elem%fcor(i,j)= 2*omega* &
+                elem%fcor(i,j)= 2*omega* &
                      (-cos(lon)*cos(lat)*sin(rangle) + sin(lat)*cos(rangle))
              else
                 elem%fcor(i,j) = 2.0D0*omega*SIN(elem%spherep(i,j)%lat)
              endif
+#ifdef HOMMEDA
+             !DA term, only nonrotated version
+             elem%fcorcosine(i,j) = 2.0D0*omega*COS(elem%spherep(i,j)%lat)
+#endif
        end do
     end do
 

components/homme/src/share/element_mod.F90

@@ -128,8 +128,10 @@ module element_mod
 
      integer(kind=long_kind)  :: gdofP(np,np)                         ! global degree of freedom (P-grid)
 
-     real (kind=real_kind)    :: fcor(np,np)                          ! Coriolis term
-
+     real (kind=real_kind)    :: fcor(np,np)                          ! Coriolis term sin
+#ifdef HOMMEDA
+     real (kind=real_kind)    :: fcorcosine(np,np)                    ! Coriolis term cos
+#endif
      type (index_t) :: idxP
      integer :: FaceNum
 

components/homme/src/share/namelist_mod.F90

@@ -819,7 +819,9 @@ subroutine readnl(par)
     call MPI_bcast(smooth_phis_nudt,1,MPIreal_t   ,par%root,par%comm,ierr)
     call MPI_bcast(initial_total_mass ,1,MPIreal_t   ,par%root,par%comm,ierr)
     call MPI_bcast(u_perturb     ,1,MPIreal_t   ,par%root,par%comm,ierr)
+
     call MPI_bcast(rotate_grid   ,1,MPIreal_t   ,par%root,par%comm,ierr)
+
     call MPI_bcast(integration,MAX_STRING_LEN,MPIChar_t ,par%root,par%comm,ierr)
 #ifndef HOMME_WITHOUT_PIOLIBRARY
     call MPI_bcast(mesh_file,MAX_FILE_LEN,MPIChar_t ,par%root,par%comm,ierr)
@@ -911,6 +913,22 @@ subroutine readnl(par)
     call MPI_bcast(use_column_solver, 1, MPIlogical_t, par%root, par%comm, ierr)
 #endif
 
+!all DA warnings together
+#ifdef HOMMEDA
+    if ( rotate_grid /= 0 ) then
+       if(par%masterproc) print *, 'DA (deep atmosphere) cannot run with rotate_grid != 0'
+       call abortmp('stopping')
+    endif
+    if ( ( rsplit == 0 ) .or. (dt_remap_factor == 0) ) then
+       if(par%masterproc) print *, 'DA (deep atmosphere) cannot run with rsplit or dt_remap = 0'
+       call abortmp('stopping')
+    endif
+    if ( theta_hydrostatic_mode ) then
+       if(par%masterproc) print *, 'DA (deep atmosphere) cannot run with theta_hydrostatic_mode=T'
+       call abortmp('stopping')
+    endif
+#endif
+
     ! should we assume Q(:,:,:,1) has water vapor:
     use_moisture = ( moisture /= "dry") 
     if (qsize<1) use_moisture = .false.  
@@ -1271,6 +1289,12 @@ subroutine readnl(par)
        write(iulog,*)"physconst: rearth = ",rearth
        write(iulog,*)"physconst: rrearth= ",rrearth
        write(iulog,*)""
+#ifdef HOMMEDA
+       write(iulog,*)"--- DEEP ATMOSPHERE ----------"
+#else
+       write(iulog,*)"--- SHALLOW ATMOSPHERE -------"
+#endif
+       write(iulog,*)""
 #endif
 
 #if !defined(CAM) && !defined(SCREAM)