Merge pull request #277 from chihchen24/cam_development_contrail

cam6_3_014: contrail parameterization added

Author: cacraigucar

cime_config/testdefs/testlist_cam.xml

@@ -1064,6 +1064,14 @@
       <option name="wallclock">00:20:00</option>
     </options>
   </test>
+  <test compset="FSD" grid="f09_f09_mg17" name="ERP_D_Ln9" testmods="cam/outfrq9s_contrail">
+    <machines>
+      <machine name="cheyenne" compiler="intel" category="prealpha"/>
+    </machines>
+    <options>
+      <option name="wallclock">00:20:00</option>
+    </options>
+  </test>
   <test compset="FSD" grid="f09_f09_mg17" name="SMS_D_Ln9" testmods="cam/outfrq9s">
     <machines>
       <machine name="cheyenne" compiler="intel" category="aux_cam"/>

cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/shell_commands

@@ -0,0 +1,2 @@
+./xmlchange ROF_NCPL=\$ATM_NCPL
+./xmlchange GLC_NCPL=\$ATM_NCPL

cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_cam

@@ -0,0 +1,8 @@
+mfilt=1,1,1,1,1,1
+ndens=1,1,1,1,1,1
+nhtfrq=9,9,9,9,9,9
+inithist='ENDOFRUN'
+aircraft_datapath ='$DIN_LOC_ROOT/atm/cam/ggas'
+aircraft_specifier              =
+ 'ac_H2O         ->   ac_H2O_filelist_monthly_2006.txt'
+ 'ac_SLANT_DIST  ->   ac_SLANT_DIST_filelist_monthly_2006.txt'

cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_clm

@@ -0,0 +1,27 @@
+!----------------------------------------------------------------------------------
+! Users should add all user specific namelist changes below in the form of 
+! namelist_var = new_namelist_value 
+!
+! Include namelist variables for drv_flds_in ONLY if -megan and/or -drydep options
+! are set in the CLM_NAMELIST_OPTS env variable.
+!
+! EXCEPTIONS: 
+! Set use_cndv           by the compset you use and the CLM_BLDNML_OPTS -dynamic_vegetation setting
+! Set use_vichydro       by the compset you use and the CLM_BLDNML_OPTS -vichydro           setting
+! Set use_cn             by the compset you use and CLM_BLDNML_OPTS -bgc  setting
+! Set use_crop           by the compset you use and CLM_BLDNML_OPTS -crop setting
+! Set spinup_state       by the CLM_BLDNML_OPTS -bgc_spinup      setting
+! Set irrigate           by the CLM_BLDNML_OPTS -irrig           setting
+! Set dtime              with L_NCPL                             option
+! Set fatmlndfrc         with LND_DOMAIN_PATH/LND_DOMAIN_FILE    options
+! Set finidat            with RUN_REFCASE/RUN_REFDATE/RUN_REFTOD options for hybrid or branch cases
+!                        (includes $inst_string for multi-ensemble cases)
+! Set glc_grid           with CISM_GRID                          option
+! Set glc_smb            with GLC_SMB                            option
+! Set maxpatch_glcmec    with GLC_NEC                            option
+! Set glc_do_dynglacier  with GLC_TWO_WAY_COUPLING               env variable
+!----------------------------------------------------------------------------------
+hist_nhtfrq = 9
+hist_mfilt  = 1
+hist_ndens  = 1
+

doc/ChangeLog

@@ -1,3 +1,69 @@
+
+===============================================================
+
+Tag name: cam6_3_014
+Originator(s): cchen, cacraig
+Date: March 23, 2021
+One-line Summary: Contrail parametrization added to CAM
+Github PR URL: https://github.com/ESCOMP/CAM/pull/277
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+Contrail Parameterization:  https://github.com/ESCOMP/CAM/issues/274
+
+Describe any changes made to build system:
+
+Describe any changes made to the namelist:
+ - "aircraft_specifier" can now take the values ac_SLANT_DIST and/or ac_TRACK_DIST to signal a contrail run
+
+List any changes to the defaults for the boundary datasets:
+
+Describe any substantial timing or memory changes:
+
+Code reviewed by: cacraig, goldy, fvitt, nusbaume
+
+List all files eliminated:
+
+List all files added and what they do:
+A       cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/shell_commands
+A       cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_cam
+A       cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_clm
+          - introduce test for contrail
+
+A       src/physics/cam/ssatcontrail.F90
+          - Actual code to process contrails
+
+List all existing files that have been modified, and describe the changes:
+M       cime_config/testdefs/testlist_cam.xml
+          - Add prealpha contrail test
+
+M       src/chemistry/utils/aircraft_emit.F90
+          - Adds logic to setup contrail run if either or both of ac_SLANT_DIST and ac_TRACK_DIST are in 
+            the aircraft_specifier list
+
+M       src/chemistry/utils/horizontal_interpolate.F90
+          - add additional logic for using flight distance and making additional calculations
+
+M       src/chemistry/utils/tracer_data.F90
+          - add additional logic if a contrail run is specified
+
+M       src/physics/cam/physpkg.F90
+          - add call to ssatcontrail_d0 (which returns immediately if not a contrail run)
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below.  All failed tests must be justified.
+
+cheyenne/intel/aux_cam: all BFB
+
+izumi/nag/aux_cam: all BFB except:
+  DAE.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall: FAIL) details:
+    FAIL DAE.f45_f45_mg37.FHS94.izumi_nag.cam-dae COMPARE_base_da
+    FAIL ERC_D_Ln9.f10_f10_mg37.QPC5.izumi_nag.cam-outfrq3s_carma TPUTCOMP Error: TPUTCOMP: Computation time increase > 25% from baseline
+     - Known failure from previous tag
+
+izumi/pgi/aux_cam: all BFB
+===============================================================
 ===============================================================
 
 Tag name: cam6_3_013

src/chemistry/utils/aircraft_emit.F90

@@ -23,6 +23,7 @@ module aircraft_emit
   public :: aircraft_emit_adv
   public :: aircraft_emit_register
   public :: aircraft_emit_readnl
+  public :: get_aircraft
 
   type :: forcing_air
      real(r8)              :: mw
@@ -40,16 +41,15 @@ module aircraft_emit
 
   type(forcing_air),allocatable :: forcings_air(:)
 
-  integer, parameter :: N_AERO = 10 
-  character(len=11)    :: aero_names(N_AERO) = (/'ac_HC      ','ac_NOX     ','ac_PMNV    ',&
-                          'ac_PMSO    ','ac_PMFO    ','ac_FUELBURN','ac_CO2     ','ac_H2O     ',&
-                          'ac_SOX     ','ac_CO      '/)
+  integer, parameter :: N_AERO = 13
+  character(len=13)    :: aero_names(N_AERO) = (/'ac_SLANT_DIST','ac_TRACK_DIST','ac_HC        ','ac_NOX       ','ac_PMNV      ',&
+                                                 'ac_PMSO      ','ac_PMFO      ','ac_FUELBURN  ','ac_CO2       ','ac_H2O       ',&
+                                                 'ac_SOX       ','ac_CO        ','ac_BC        '/)
 
   real(r8), parameter :: molmass(N_AERO) = 1._r8
 
-  logical :: advective_tracer(N_AERO) = (/.false., .false., .false., .false., .false., &
-                                          .false., .false., .false., .false.,.false./)
-  character(len=3) :: mixtype(N_AERO) = (/'wet','wet','wet','wet','wet','wet','wet','wet','wet','wet'/)
+  logical :: advective_tracer(N_AERO) = .false.
+  character(len=3) :: mixtype(N_AERO) = 'wet'
 
   real(r8) :: cptmp = 666.0_r8
   real(r8) :: qmin = 0.0_r8
@@ -67,9 +67,26 @@ module aircraft_emit
   integer :: aircraft_cnt = 0
   character(len=16) :: spc_name_list(N_AERO)
   character(len=256) :: spc_flist(N_AERO),spc_fname(N_AERO)
+  logical :: dist(N_AERO)
 
 contains
 
+  subroutine get_aircraft(cnt, spc_name_list_out)
+  integer, intent(out) :: cnt 
+  character(len=16), optional, intent(out) :: spc_name_list_out(N_AERO)
+  integer :: i
+
+  spc_name_list_out = ''
+
+  cnt = aircraft_cnt
+  if( cnt>0 ) then
+    do i=1,cnt
+       spc_name_list_out(i) = spc_name_list(i)
+    end do
+  end if
+
+  end subroutine get_aircraft
+ 
   subroutine aircraft_emit_register()
 
 !------------------------------------------------------------------
@@ -90,6 +107,8 @@ subroutine aircraft_emit_register()
     !------------------------------------------------------------------
     ! Return if air_specifier is blank (no aircraft data to process)
     !------------------------------------------------------------------
+    dist(:) = .false.
+    aircraft_cnt = 0
     if (air_specifier(1) == "") return
 
 ! count aircraft emission species used in the simulation
@@ -108,6 +127,8 @@ subroutine aircraft_emit_register()
 	 call endrun('aircraft_emit_register: '//trim(spc_name)//' is not in the aircraft emission dataset')
         endif
 
+        if (trim(spc_name) == 'ac_SLANT_DIST'.or. trim(spc_name) == 'ac_TRACK_DIST') dist(n) = .true.
+
         aircraft_cnt = aircraft_cnt + 1
         call pbuf_add_field(aero_names(mm),'physpkg',dtype_r8,(/pcols,pver/),idx)
 
@@ -189,6 +210,7 @@ subroutine aircraft_emit_init()
           forcings_air(m)%file%cyclical_list    = .true.  ! Aircraft data cycles over the filename list
           forcings_air(m)%file%weight_by_lat     = .true.  ! Aircraft data -  interpolated with latitude weighting
           forcings_air(m)%file%conserve_column = .true. ! Aircraft data - vertically interpolated to conserve the total column
+          forcings_air(m)%file%dist = dist(m)
           forcings_air(m)%species          = spc_name
           forcings_air(m)%sectors          = spc_name ! Only one species per file for aircraft data
           forcings_air(m)%nsectors         = 1
@@ -314,6 +336,8 @@ subroutine aircraft_emit_adv( state, pbuf2d)
               caseid = 4
           case ('kg m-2 s-1')
               caseid = 5
+          case ('m/sec' ) 
+              caseid = 6
           case default
              print*, 'aircraft_emit_adv: units = ',trim(forcings_air(m)%fields(i)%units) ,' are not recognized'
              call endrun('aircraft_emit_adv: units are not recognized')
@@ -339,6 +363,8 @@ subroutine aircraft_emit_adv( state, pbuf2d)
                 to_mmr(:ncol,:) = 1.0_r8
              elseif (caseid == 5) then
                 to_mmr(:ncol,:) = 1.0_r8
+             elseif (caseid == 6) then
+                to_mmr(:ncol,:) = 1.0_r8
              else
                 to_mmr(:ncol,:) = molmass(ind)/mwdry
              endif

src/chemistry/utils/horizontal_interpolate.F90

@@ -17,7 +17,7 @@ module horizontal_interpolate
   public :: xy_interp_init, xy_interp
 
 contains
-  subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y)
+  subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight_distance)
 !------------------------------------------------------------------------------------------------------------
 ! This program computes weighting functions to map a variable of (im1,jm1) resolution to (im2,jm2) resolution
 ! weight_x(im2,im1) is the weighting function for zonal interpolation
@@ -28,6 +28,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y)
 !------------------------------------------------------------------------------------------------------------
   implicit none
   integer,  intent(in)  :: im1, jm1, im2, jm2
+  logical,  intent(in)  :: use_flight_distance   !.true. = flight distance, .false. =  all mixing ratios
   real(r8), intent(in)  :: lon0(im1), lat0(jm1)
   real(r8), intent(out) :: weight_x(im2,im1), weight_y(jm2,jm1)
 
@@ -110,28 +111,40 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y)
 ! check if there is any overlap between the source grid and the target grid
 ! if no overlap, then weighting is zero
 ! there are three scenarios overlaps can take place 
-   if( (x1_west.ge.x2_west).and.(x1_east.le.x2_east) ) then
+   if( (x1_west>=x2_west).and.(x1_east<=x2_east) ) then
 ! case 1: 
 !                x1_west             x1_east
 !                  |-------------------|
 !            |---------------------------------|
 !          x2_west                           x2_east
+    if(use_flight_distance) then
+     weight_x(i2,i1) = 1.0_r8
+    else
      weight_x(i2,i1) =  (x1_east-x1_west)/(x2_east-x2_west)
-   elseif ( (x1_west.ge.x2_west).and.(x1_west.lt.x2_east) ) then
+    endif
+   elseif ( (x1_west>=x2_west).and.(x1_west<x2_east) ) then
 ! case 2: 
 !                x1_west                          x1_east
 !                  |--------------------------------|
 !            |---------------------------------|
 !          x2_west                           x2_east
+   if(use_flight_distance) then
+     weight_x(i2,i1) = (x2_east-x1_west)/(x1_east-x1_west)
+   else
      weight_x(i2,i1) = (x2_east-x1_west)/(x2_east-x2_west)
-   elseif ( (x1_east>x2_west).and.(x1_east.le.x2_east) ) then
+   endif
+   elseif ( (x1_east>x2_west).and.(x1_east<=x2_east) ) then
 ! case 3: 
 !       x1_west                          x1_east
 !         |--------------------------------|
 !                |---------------------------------|
 !              x2_west                           x2_east
+   if(use_flight_distance) then
+     weight_x(i2,i1) = (x1_east-x2_west)/(x1_east-x1_west)
+   else
      weight_x(i2,i1) = (x1_east-x2_west)/(x2_east-x2_west)
-   elseif ( (x1_east.gt.x2_east).and.(x1_west.lt.x2_west) ) then
+   endif
+   elseif ( (x1_east>x2_east).and.(x1_west<x2_west) ) then
 ! case 4: 
 !       x1_west                          x1_east
 !         |--------------------------------|
@@ -145,22 +158,30 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y)
 
 
 ! consider end points
-      if(slon1(im1+1).gt.slon2(im2+1)) then
+      if(slon1(im1+1)>slon2(im2+1)) then
 ! case 1:
 !           slon1(im1)                slon1(im1+1) <--- end point
 !              |-------------------------|
 !           |----------------|......................|
 !        slon2(im2)         slon2(im2+1)        slon2(2)  (note: slon2(im2+1) = slon2(1))
-     	weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon2(2)-slon2(1))
+         if(use_flight_distance) then
+            weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon1(im1+1)-slon1(im1))
+         else
+            weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon2(2)-slon2(1))
+         endif
       endif	
 
-      if(slon1(im1+1).lt.slon2(im2+1)) then
+      if(slon1(im1+1)<slon2(im2+1)) then
 ! case 1:
 !           slon1(im1)                slon1(im1+1)                  slon1(2)    (note: slon1(im1+1) = slon1(1))
 !              |-------------------------|.............................|
 !                   |-------------------------------|
 !               slon2(im2)                        slon2(im2+1) <--- end point
-        weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon2(2)-slon2(1)) 
+         if(use_flight_distance) then
+            weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon1(2)-slon1(1))
+         else
+            weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon2(2)-slon2(1)) 
+         endif
       endif
 
 
@@ -181,34 +202,50 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y)
 ! there are three scenarios overlaps can take place 
 ! note: there is Guassian weight to consider in the meridional direction!
 
-            if( (y1_south.ge.y2_south).and.(y1_north.le.y2_north) ) then
+            if( (y1_south>=y2_south).and.(y1_north<=y2_north) ) then
 ! case 1: 
 !                y1_south             y1_north
 !                  |-------------------|
 !            |---------------------------------|
 !          y2_south                           y2_north
-                 weight_y(j2,j1) =  gw1(j1)/gw2(j2)
-            elseif ( (y1_south.ge.y2_south).and.(y1_south.lt.y2_north) ) then
+                if(use_flight_distance) then
+                   weight_y(j2,j1) =  1.0_r8
+                else
+                   weight_y(j2,j1) =  gw1(j1)/gw2(j2)
+                endif
+            elseif ( (y1_south>=y2_south).and.(y1_south<y2_north) ) then
 ! case 2: 
 !                y1_south                          y1_north
 !                  |--------------------------------|
 !            |---------------------------------|
 !          y2_south                           y2_north
-                 weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2)
-            elseif ( (y1_north.gt.y2_south).and.(y1_north.le.y2_north) ) then
+                if(use_flight_distance) then
+                   weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south)
+                else
+                   weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2)
+                endif
+            elseif ( (y1_north>y2_south).and.(y1_north<=y2_north) ) then
 ! case 3: 
 !       y1_south                          y1_north
 !         |--------------------------------|
 !                |---------------------------------|
 !              y2_south                           y2_north
-                 weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2)
-            elseif ( (y1_north.gt.y2_north).and.(y1_south.lt.y2_south) ) then
+                if(use_flight_distance) then
+                   weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)
+                else
+                   weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2)
+                endif
+            elseif ( (y1_north>y2_north).and.(y1_south<y2_south) ) then
 ! case 4: 
 !       y1_south                          y1_north
 !         |--------------------------------|
 !                |---------------------|
 !              y2_south             y2_north
-                 weight_y(j2,j1) = gw1(j1)/gw2(j2)
+                if(use_flight_distance) then
+                   weight_y(j2,j1) = 1._r8
+                else
+                   weight_y(j2,j1) = gw1(j1)/gw2(j2)
+                endif
             endif
 
           enddo