Merge branch 'halldm2000/homme/dcmip12_test1.x' PR #1100

[BFB]
Added DCMIP tracer transport tests 1-1, 1-2, and 1-3
moved prescribed_velocity code into its own subroutine in prim_advance
verified that all tests 1.x work on Edison
verified the HOMME stand-alone tests pass (12/12)
verified ACME Aquaplanet test builds, runs, and looks normal

components/homme/README.cmake

@@ -1,6 +1,6 @@
 03/2013 BFJ
 02/2016 MT
-
+10/2016 DMH
 
 Please see the HOMME wiki for information on how to build HOMME using the CMake build system.
 https://wiki.ucar.edu/display/homme/The+HOMME+CMake+build+and+testing+system
@@ -29,10 +29,13 @@ $WDIR/test_execs/swtcB     directory for swtcB executable
 $WDIR/utils         cprnc utility, and PIO and timing libraries
 $WDIR/tests     directory containing all the HOMME regression tests
 
-
 HOMME has a large regression test suite.  For instructions on running and adding
 new tests, see homme/test/reg_test/README
 
+DCMIP tests provide a standard means for testing and comparing the ACME HOMME dycore with other dycores
+both hydrostatic and nonhydrostatic. They have been placed in their own dcmip_test directory for now.
+To run a DCMIP tests, navigate to the appropriate directory and type make install to install test scripts and namelists.
+
 The CMAKE code could use some cleanup. 
 - user configured variables should not need to be prefixed by the exectuable name
   (i.e. -DNP=4, instead of -DPREQX_NP=4).  This will make the cmake code a lot simpler
@@ -41,9 +44,6 @@ The CMAKE code could use some cleanup.
   with cmake's tree-like directory approach, the tests should be associated with their
 
 
-
-
-
 ************************************************************************************************
 
 ***OBSOLETE***

components/homme/dcmip_tests/dcmip2012_test1.1_3d_deformational_flow/CMakeLists.txt

@@ -0,0 +1,14 @@
+file(GLOB NCL_SCRIPTS    "*.ncl")                                       # get ncl-scripts
+file(GLOB SHELL_SCRIPTS  "*.sh")                                        # get shell-scripts
+file(GLOB NAMELISTS      "*.nl")                                        # get namelists
+file(GLOB PYTHON_SCRIPTS "*.py")                                        # get python scripts
+file(GLOB M_FILES "*.m")                                                # get matlab scripts
+
+# Copy test files to build
+install(PROGRAMS ${NCL_SCRIPTS} ${SHELL_SCRIPTS} ${NAMELISTS} ${PYTHON_SCRIPTS} ${M_FILES} DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+
+# Create a symbolic link to the vertical coordinate directory
+execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_SOURCE_DIR}/test/vcoord ${CMAKE_CURRENT_BINARY_DIR}/vcoord )
+
+# Add subdirectory for results
+add_subdirectory(movies)

components/homme/dcmip_tests/dcmip2012_test1.1_3d_deformational_flow/README.txt

@@ -0,0 +1,20 @@
+https://www.earthsystemcog.org/projects/dcmip-2012/cam-se_test_11_new
+
+Passive 3d advection with prescribed winds.
+
+Velocities reversed at the midpoint, returning tracers to initial position.
+Wind velocities prescribed. Dynamic updates of v, T, p disabled. 
+
+Test 1-1: 3D Deformational Flow
+
+Five tracer fields transported simultaneously:
+
+1. q1 = two cosine bells
+2. q2 = inverse squared cosine bells
+3. q3 = a slotted ellipse
+4. q4 = 1 - 0.3(q1+q2+q3)
+5. q5 = 1
+
+test resolution: 
+1 degree horizontal, 60 uniformly-spaced vertical levels, duration 12 days
+

components/homme/dcmip_tests/dcmip2012_test1.1_3d_deformational_flow/build-preqx.sh

@@ -0,0 +1,6 @@
+#!/bin/bash
+cwd=`pwd`
+cd ../..
+  echo "make preqx-nlev60-interp"
+  make -j4 preqx-nlev60-interp
+cd $cwd

components/homme/dcmip_tests/dcmip2012_test1.1_3d_deformational_flow/jobscript-darwin-preqx.sh

@@ -0,0 +1,8 @@
+#!/bin/bash
+#
+#   Jobscript for launching dcmip2012 test1-1 on a mac running Darwin
+#
+# usage: ./jobscript-...
+
+EXEC=../../test_execs/preqx-nlev60-interp/preqx-nlev60-interp     # set name of executable
+openmpiexec -n 6 $EXEC < ./namelist-lowres.nl       # launch simulation

components/homme/dcmip_tests/dcmip2012_test1.1_3d_deformational_flow/jobscript-edison-preqx.sh

@@ -0,0 +1,16 @@
+#!/bin/bash
+#
+#   Jobscript for launching dcmip2012 test1-1 on the NERSC Edison machine
+#
+# usage: sbatch jobscript-...
+
+#SBATCH -J dcmip1-1           # job name
+#SBATCH -o out_dcmip1-1.o%j   # output and error file name (%j expands to jobID)
+#SBATCH -n 768                # total number of mpi tasks requested
+#SBATCH -p debug              # queue (partition) -- normal, development, etc.
+#SBATCH -t 00:15:00           # run time (hh:mm:ss)
+#SBATCH -A acme               # charge hours to account 1
+
+EXEC=../../test_execs/preqx-nlev60-interp/preqx-nlev60-interp           # set name of executable
+srun -n 768 $EXEC < ./namelist-default.nl                               # launch simulation
+

components/homme/dcmip_tests/dcmip2012_test1.1_3d_deformational_flow/namelist-default.nl

@@ -0,0 +1,50 @@
+!
+! namelist for dcmip2012 test1-1: 3d deformational flow
+!_______________________________________________________________________
+&ctl_nl
+  nthreads          = 1
+  partmethod        = 4                         ! mesh parition method: 4 = space filling curve
+  topology          = "cube"                    ! mesh type: cubed sphere
+  test_case         = "dcmip2012_test1_1"       ! test identifier
+  ne                = 30                        ! number of elements per cube face
+  qsize             = 4                         ! num tracer fields
+  ndays             = 12                         ! num simulation days: 0 = use nmax steps
+  statefreq         = 200                       ! number of steps between screen dumps
+  restartfreq       = -1                        ! don't write restart files if < 0
+  runtype           = 0                         ! 0 = new run
+  tstep             = 33                        ! largest timestep in seconds
+  integration       = 'explicit'                ! explicit time integration
+  tstep_type        = 1                         ! 1 => default method
+  smooth            = 0                         ! timestep smooting
+  nu                = 1.585e13                  ! hyperviscosity minimum L2 error at 10^13.2
+  nu_s              = 1.585e13
+  hypervis_order    = 2                         ! 2 = hyperviscosity
+  hypervis_subcycle = 1                         ! 1 = no hyperviz subcycling
+  prescribed_wind   = 1
+/
+&filter_nl/
+&solver_nl
+  precon_method     = "identity"
+  maxits            = 50
+  tol               = 1.e-7
+/
+&vert_nl
+  vform             = "ccm"                     ! vertical coordinate type "ccm"=hybrid pressure/terrain
+  vanalytic         = 1                         ! set vcoords in initialization routine
+  vtop              = 0.2549944                 ! vertical coordinate at top of atm 254.9 hPa (12km)
+/
+&analysis_nl
+  output_dir        = "./movies/"               ! destination dir for netcdf file
+  output_timeunits  = 1,                        ! 1=days, 2=hours, 0=timesteps
+  output_frequency  = 1,                        ! output every 4 hours
+  output_varnames1  ='u','geo','Q','Q2','Q3','Q4' ! variables to write to file
+  interp_type       = 0                         ! 0=native grid, 1=bilinear
+  output_type       ='netcdf'                   ! netcdf or pnetcdf
+  num_io_procs      = 16         
+  interp_nlat       = 256
+  interp_nlon       = 512
+/
+&prof_inparm
+  profile_outpe_num   = 100
+  profile_single_file	= .true.
+/

components/homme/dcmip_tests/dcmip2012_test1.1_3d_deformational_flow/namelist-lowres.nl

@@ -0,0 +1,50 @@
+!
+! namelist for dcmip2012 test1-1: 3d deformational flow
+!_______________________________________________________________________
+&ctl_nl
+  nthreads          = 1
+  partmethod        = 4                         ! mesh parition method: 4 = space filling curve
+  topology          = "cube"                    ! mesh type: cubed sphere
+  test_case         = "dcmip2012_test1_1"       ! test identifier
+  ne                = 7                         ! number of elements per cube face
+  qsize             = 4                         ! num tracer fields
+  ndays             = 12                        ! num simulation days: 0 = use nmax steps
+  statefreq         = 10 !50                        ! number of steps between screen dumps
+  restartfreq       = -1                        ! don't write restart files if < 0
+  runtype           = 0                         ! 0 = new run
+  tstep             = 100 !33                        ! largest timestep in seconds
+  integration       = 'explicit'                ! explicit time integration
+  tstep_type        = 1                         ! 1 => default method
+  smooth            = 0                         ! timestep smooting
+  nu                = 1.585e13                  ! hyperviscosity minimum L2 error at 10^13.2
+  nu_s              = 1.585e13
+  hypervis_order    = 2                         ! 2 = hyperviscosity
+  hypervis_subcycle = 1                         ! 1 = no hyperviz subcycling
+  prescribed_wind   = 1
+/
+&filter_nl/
+&solver_nl
+  precon_method     = "identity"
+  maxits            = 50
+  tol               = 1.e-7
+/
+&vert_nl
+  vform             = "ccm"                     ! vertical coordinate type "ccm"=hybrid pressure/terrain
+  vanalytic         = 1                         ! set vcoords in initialization routine
+  vtop              = 0.2549944                 ! vertical coordinate at top of atm 254.9 hPa (12km)
+/
+&analysis_nl
+  output_dir        = "./movies/"               ! destination dir for netcdf file
+  output_timeunits  = 2,                        ! 1=days, 2=hours, 0=timesteps
+  output_frequency  = 4,                        ! output every 4 hours
+  output_varnames1  ='u','geo','Q','Q2','Q3','Q4' ! variables to write to file
+  interp_type       = 0                         ! 0=native grid, 1=bilinear
+  output_type       ='netcdf'                   ! netcdf or pnetcdf
+  num_io_procs      = 16         
+  interp_nlat       = 256
+  interp_nlon       = 512
+/
+&prof_inparm
+  profile_outpe_num   = 100
+  profile_single_file	= .true.
+/

components/homme/dcmip_tests/dcmip2012_test1.1_3d_deformational_flow/plot_q_lat_lon.ncl

@@ -0,0 +1,186 @@
+;================================================;
+;  plot-q-lat-lon
+;  plot lattitude-longitude cross section of tracer q
+;  at time_index t
+;================================================;
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"   
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   
+; ================================================;
+
+begin
+
+  f    = addfile("movies/dcmip2012_test1_11.nc","r")
+  lat  = f->lat
+  lon  = f->lon
+  lev  = f->lev
+  time = f->time
+  nlat = getfilevardimsizes(f, "lat" )
+  nlon = getfilevardimsizes(f, "lon" )
+  nlev = getfilevardimsizes(f, "lev" )
+  nt   = dimsizes(time)
+
+  ; get time_index t from command line
+
+  t_index = nt-1
+  if(isvar("t")) then
+    t_index = t
+  end if
+
+  ; get tracer_id q from command line
+
+  tracer_id=1
+  if(isvar("q")) then
+    tracer_id = q         ; use q from command line if possible
+  end if
+
+  ; get zoom factor from command line
+
+  zoom_factor = 1.0
+  if(isvar("zoom"))
+    zoom_factor = zoom
+  end if
+
+  ; find level closest to 4900 m
+
+  g = 9.80616d0; Gravity (m s^2)
+  heights = f->geo(0,:,0,0)/g
+;  level = (nlev-1) - closest_val(4900,heights(::-1))
+  level = (nlev-1) - closest_val(5100,heights(::-1))
+
+  ; use level index from command line if specified
+  if(isvar("l"))
+    level = l
+  end if
+
+  print("height = "+heights(level))
+  z = round( heights(level), 3)
+
+  print ("tracer_id:    q    = " + tracer_id)
+  print ("t_index:      t    = " + t_index)
+  print ("zoom_factor:  zoom = " + zoom_factor)
+  print ("level:        l    = " + level)
+  print ("eta:               = " + lev(level) )
+  print ("z:                 = " + z )
+
+  ; load chosen field
+
+  if(tracer_id .eq. 1) then
+    q_i = f->Q (t_index,level,:,:)
+
+  else if (tracer_id .eq. 2) then
+    q_i = f->Q2(t_index,level,:,:)
+
+  else if (tracer_id .eq. 3) then
+    q_i = f->Q3(t_index,level,:,:)
+
+  else if (tracer_id .eq. 4) then
+    q_i = f->Q4(t_index,level,:,:)
+
+  end if
+  end if
+  end if
+  end if
+
+  ; delete units and long names
+
+  q_i@units = " "
+  q_i@long_name = " "
+
+  ; display min and max values
+  minval = min(q_i)
+  maxval = max(q_i)
+  print("minval="+minval+" maxval="+maxval)
+
+  ; round to two decimal places?
+  ;q_i = decimalPlaces(q_i,2,True)
+  ;printMinMax(q_i,0)
+
+  ; Setup the plot
+
+  res1                      = True
+  res1@gsnDraw              = False        		; panel plot
+  res1@gsnFrame             = False           ; don't draw yet
+  res1@gsnCenterString      = ""
+  res1@vpWidthF             = 2
+  res1@vpHeightF            = 1
+  res1@sfXArray             = lon           	; uses lon as plot x-axis
+  res1@sfYArray             = lat         		; uses lat for y axis
+  res1@tiYAxisString        = "latitiude"     ; y-axis title
+  res1@tiXAxisString        = "longitude"      ; x-axis title
+  res1@lbLabelBarOn         = False           ; label bar
+
+  res1@cnLinesOn            = False ;True
+  res1@cnFillOn             = True
+  res1@cnInfoLabelOn        = False           ; don't give each
+  res1@cnLineLabelsOn       = False           ; panel its own
+  res1@cnMonoLineColor      = True           ; single color lines
+  res1@cnLineColor          = "black"         ; black lines
+  res1@cnLineDrawOrder      = "PreDraw"
+  res1@cnLineThicknessF     = 0.1
+  res1@cnFillOpacityF       = 1.0
+  res1@cnLevelSelectionMode = "ManualLevels"
+  res1@cnMaxLevelValF       = 1.00            ; max contour color label
+  res1@cnMinLevelValF       = 0.00            ; min contour color label
+  res1@cnLevelSpacingF      = 0.10            ; contour color spacing
+
+  ; text properties
+
+  res1@tmXBMajorLengthF     = 0.03
+  res1@txFont               = 22              ; helvetica-bold
+  res1@tiMainFont           = 22              ; helvetica-bold
+  res1@tiMainFontHeightF    = 0.045            ; set font sizes
+  res1@tiXAxisFontHeightF   = 0.05
+  res1@tiYAxisFontHeightF   = 0.05
+  res1@tmXBLabelFontHeightF = 0.04
+  res1@tiXAxisOffsetYF      = 0.010           ; shift x label
+  res1@tiYAxisOffsetXF      = 0.010           ; shift y label
+	res1@tiMainString         = "DCMIP 1-1  q"+tracer_id+"  t="+sprintf("%3.1f", time(t_index))+" days"+ "  z="+z+"m  range=("+sprintf("%5.4f", minval)+","+sprintf("%5.4f", maxval)+")"
+
+  ; zoom in on the center of the plot
+
+  res1@trXMinF =  180.0 - 180.0/zoom_factor
+  res1@trXMaxF =  180.0 + 180.0/zoom_factor
+  res1@trYMinF =    0.0 -  90.0/zoom_factor
+  res1@trYMaxF =    0.0 +  90.0/zoom_factor
+
+  ; panel properties
+
+  pres                      = True            ; Panel resource
+  pres@gsnMaximize          = True            ; maximize panel size
+  pres@gsnFrame             = False           ; don't advance yet
+  pres@gsnPaperOrientation  = "portrait"      ; paper orientation
+  pres@gsnPanelLabelBar     = True            ; add common label bar
+  pres@lbLabelFontHeightF   = 0.020           ; set labelbar font height
+
+  ; set output type and name
+
+  wks_type                  = "pdf"
+  wks_type@wkPaperHeightF   = 8
+  wks_type@wkPaperWidthF    = 11
+  filename ="dcmip2012_test1-1_q"+tracer_id+"_t"+t_index
+	wks  = gsn_open_wks(wks_type,filename)
+  print("writing "+filename)
+
+  ; use rainbow colormap with a few of color substitutions
+
+  gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")      ; choose default colormap
+  colors       = gsn_retrieve_colormap(wks)   ; retrieve color map for editing
+  colors(25 ,:)= (/ 1, 1, 1/)                 ; white     for tracer-free region
+  colors(2  ,:)= (/ 0.95, 0.95, 0.95 /)          ; gray for undershoots
+  colors(255,:)= (/ 0.4, 0.0, 0.0 /)          ; brick-red for overshoots
+  gsn_define_colormap(wks,colors)
+
+  ; plot countours and label bar
+
+  plot    = new (1, graphic)               			; define plot - need 4 panels
+  plot(0) = gsn_csm_contour(wks,q_i(:,:),res1)		; plot Q1
+  gsn_panel(wks,plot,(/1,1/),pres)			; 2x2 plot
+  frame(wks)
+
+end
+
+
+
+