diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml
index 0abda7822b..8ada83654b 100644
--- a/cime_config/testdefs/testlist_cam.xml
+++ b/cime_config/testdefs/testlist_cam.xml
@@ -1064,6 +1064,14 @@
+
+
+
+
+
+
+
+
diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/shell_commands
new file mode 100644
index 0000000000..eb40ad83e0
--- /dev/null
+++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/shell_commands
@@ -0,0 +1,2 @@
+./xmlchange ROF_NCPL=\$ATM_NCPL
+./xmlchange GLC_NCPL=\$ATM_NCPL
diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_cam
new file mode 100644
index 0000000000..80bb5f2bf2
--- /dev/null
+++ b/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'
diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_clm b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_clm
new file mode 100644
index 0000000000..0d83b5367b
--- /dev/null
+++ b/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
+
diff --git a/doc/ChangeLog b/doc/ChangeLog
index 8536cdd709..0a7d413f5f 100644
--- a/doc/ChangeLog
+++ b/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
diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90
index 1403a4e264..f2a242aa63 100644
--- a/src/chemistry/utils/aircraft_emit.F90
+++ b/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
diff --git a/src/chemistry/utils/horizontal_interpolate.F90 b/src/chemistry/utils/horizontal_interpolate.F90
index e19ee36aa3..6a7ae7e4f8 100644
--- a/src/chemistry/utils/horizontal_interpolate.F90
+++ b/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_westx2_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_westslon2(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)=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_southy2_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 null(), weight_y => null()
integer, pointer, dimension(:) :: count_x => null(), count_y => null()
integer, pointer, dimension(:,:) :: index_x => null(), index_y => null()
+
+ real(r8), pointer, dimension(:,:) :: weight0_x=>null(), weight0_y=>null()
+ integer, pointer, dimension(:) :: count0_x=>null(), count0_y=>null()
+ integer, pointer, dimension(:,:) :: index0_x=>null(), index0_y=>null()
+ logical :: dist
+
real(r8) :: p0
type(var_desc_t) :: ps_id
logical, allocatable, dimension(:) :: in_pbuf
@@ -162,9 +168,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
use phys_grid, only : get_rlat_all_p, get_rlon_all_p, get_ncols_p
use dycore, only : dycore_is
use horizontal_interpolate, only : xy_interp_init
-#if ( defined SPMD )
- use mpishorthand, only: mpicom, mpir8, mpiint
-#endif
+ use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_real8, mpi_integer
implicit none
@@ -180,6 +184,8 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
integer, intent(in) :: data_fixed_tod
character(len=*), intent(in) :: data_type
+ character(len=*), parameter :: sub = 'trcdata_init'
+
integer :: f, mxnflds, astat
integer :: str_yr, str_mon, str_day
integer :: lon_dimid, lat_dimid, lev_dimid, tim_dimid, old_dimid
@@ -197,6 +203,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
real(r8):: rlats(pcols), rlons(pcols)
integer :: lchnk, ncol, icol, i,j
logical :: found
+ integer :: aircraft_cnt
call specify_fields( specifier, flds )
@@ -215,6 +222,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
file%fill_in_months = .false.
file%cyclical = .false.
file%cyclical_list = .false.
+ file%dist = .false.
select case ( data_type )
case( 'FIXED' )
@@ -588,8 +596,8 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
call get_horiz_grid_d(plat, clat_d_out=phi)
call get_horiz_grid_d(plon, clon_d_out=lam)
- allocate(lon_global_grid_ndx(pcols,begchunk:endchunk))
- allocate(lat_global_grid_ndx(pcols,begchunk:endchunk))
+ if(.not.allocated(lon_global_grid_ndx)) allocate(lon_global_grid_ndx(pcols,begchunk:endchunk))
+ if(.not.allocated(lat_global_grid_ndx)) allocate(lat_global_grid_ndx(pcols,begchunk:endchunk))
lon_global_grid_ndx=-huge(1)
lat_global_grid_ndx=-huge(1)
@@ -617,12 +625,36 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
deallocate(phi,lam)
! weight_x & weight_y are weighting function for x & y interpolation
- allocate(file%weight_x(plon,file%nlon))
- allocate(file%weight_y(plat,file%nlat))
- allocate(file%count_x(plon))
- allocate(file%count_y(plat))
- allocate(file%index_x(plon,file%nlon))
- allocate(file%index_y(plat,file%nlat))
+ allocate(file%weight_x(plon,file%nlon), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%weight_x allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate weight_x array')
+ end if
+ allocate(file%weight_y(plat,file%nlat), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%weight_y allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate weight_y array')
+ end if
+ allocate(file%count_x(plon), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%count_x allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate count_x array')
+ end if
+ allocate(file%count_y(plat), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%count_y allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate count_y array')
+ end if
+ allocate(file%index_x(plon,file%nlon), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%index_x allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate index_x array')
+ end if
+ allocate(file%index_y(plat,file%nlat), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%index_y allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate index_y array')
+ end if
file%weight_x(:,:) = 0.0_r8
file%weight_y(:,:) = 0.0_r8
file%count_x(:) = 0
@@ -630,14 +662,54 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
file%index_x(:,:) = 0
file%index_y(:,:) = 0
+ if( file%dist ) then
+ allocate(file%weight0_x(plon,file%nlon), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%weight0_x allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate weight0_x array')
+ end if
+ allocate(file%weight0_y(plat,file%nlat), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%weight0_y allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate weight0_y array')
+ end if
+ allocate(file%count0_x(plon), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%count0_x allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate count0_x array')
+ end if
+ allocate(file%count0_y(plat), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%count0_y allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate count0_y array')
+ end if
+ allocate(file%index0_x(plon,file%nlon), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%index0_x allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate index0_x array')
+ end if
+ allocate(file%index0_y(plat,file%nlat), stat=astat)
+ if( astat /= 0 ) then
+ write(iulog,*) 'trcdata_init: file%index0_y allocation error = ',astat
+ call endrun('trcdata_init: failed to allocate index0_y array')
+ end if
+ file%weight0_x(:,:) = 0.0_r8
+ file%weight0_y(:,:) = 0.0_r8
+ file%count0_x(:) = 0
+ file%count0_y(:) = 0
+ file%index0_x(:,:) = 0
+ file%index0_y(:,:) = 0
+ endif
+
if(masterproc) then
! compute weighting
- call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats,plon,plat,file%weight_x,file%weight_y)
+ call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats, &
+ plon,plat,file%weight_x,file%weight_y,file%dist)
do i2=1,plon
file%count_x(i2) = 0
do i1=1,file%nlon
- if(file%weight_x(i2,i1).gt.0.0_r8 ) then
+ if(file%weight_x(i2,i1)>0.0_r8 ) then
file%count_x(i2) = file%count_x(i2) + 1
file%index_x(i2,file%count_x(i2)) = i1
endif
@@ -647,22 +719,65 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, &
do j2=1,plat
file%count_y(j2) = 0
do j1=1,file%nlat
- if(file%weight_y(j2,j1).gt.0.0_r8 ) then
+ if(file%weight_y(j2,j1)>0.0_r8 ) then
file%count_y(j2) = file%count_y(j2) + 1
file%index_y(j2,file%count_y(j2)) = j1
endif
enddo
enddo
- endif
-#if ( defined SPMD)
- call mpibcast(file%weight_x, plon*file%nlon, mpir8 , 0, mpicom)
- call mpibcast(file%weight_y, plat*file%nlat, mpir8 , 0, mpicom)
- call mpibcast(file%count_x, plon, mpiint , 0, mpicom)
- call mpibcast(file%count_y, plat, mpiint , 0, mpicom)
- call mpibcast(file%index_x, plon*file%nlon, mpiint , 0, mpicom)
- call mpibcast(file%index_y, plat*file%nlat, mpiint , 0, mpicom)
-#endif
+ if( file%dist ) then
+ call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats,&
+ plon,plat,file%weight0_x,file%weight0_y,file%dist)
+
+ do i2=1,plon
+ file%count0_x(i2) = 0
+ do i1=1,file%nlon
+ if(file%weight0_x(i2,i1)>0.0_r8 ) then
+ file%count0_x(i2) = file%count0_x(i2) + 1
+ file%index0_x(i2,file%count0_x(i2)) = i1
+ endif
+ enddo
+ enddo
+
+ do j2=1,plat
+ file%count0_y(j2) = 0
+ do j1=1,file%nlat
+ if(file%weight0_y(j2,j1)>0.0_r8 ) then
+ file%count0_y(j2) = file%count0_y(j2) + 1
+ file%index0_y(j2,file%count0_y(j2)) = j1
+ endif
+ enddo
+ enddo
+ endif
+ endif
+
+ call mpi_bcast(file%weight_x, plon*file%nlon, mpi_real8 , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight_x")
+ call mpi_bcast(file%weight_y, plat*file%nlat, mpi_real8 , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight_y")
+ call mpi_bcast(file%count_x, plon, mpi_integer , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count_x")
+ call mpi_bcast(file%count_y, plat, mpi_integer , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count_y")
+ call mpi_bcast(file%index_x, plon*file%nlon, mpi_integer , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index_x")
+ call mpi_bcast(file%index_y, plat*file%nlat, mpi_integer , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index_y")
+ if( file%dist ) then
+ call mpi_bcast(file%weight0_x, plon*file%nlon, mpi_real8 , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight0_x")
+ call mpi_bcast(file%weight0_y, plat*file%nlat, mpi_real8 , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight0_y")
+ call mpi_bcast(file%count0_x, plon, mpi_integer , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count0_x")
+ call mpi_bcast(file%count0_y, plon, mpi_integer , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count0_y")
+ call mpi_bcast(file%index0_x, plon*file%nlon, mpi_integer , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index0_x")
+ call mpi_bcast(file%index0_y, plat*file%nlat, mpi_integer , mstrid, mpicom,ierr)
+ if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index0_y")
+ endif
endif
end subroutine trcdata_init
@@ -1087,7 +1202,7 @@ subroutine find_times( itms, fids, time, file, datatimem, datatimep, times_found
datatimep = all_data_times(np1) !+ file%offset_time
! When stepTime, datatimep may not equal the time (as only datatimem is used)
! Should not break other runs?
- if ( (time .ge. datatimem) .and. (time .lt. datatimep) ) then
+ if ( (time >= datatimem) .and. (time < datatimep) ) then
times_found = .true.
exit find_times_loop
endif
@@ -1626,16 +1741,27 @@ subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order)
if(file%weight_by_lat) then
call t_startf('xy_interp')
-
- do c = begchunk,endchunk
+ if( file%dist ) then
+ do c = begchunk,endchunk
ncols = get_ncols_p(c)
lons(:ncols) = lon_global_grid_ndx(:ncols,c)
lats(:ncols) = lat_global_grid_ndx(:ncols,c)
- call xy_interp(file%nlon,file%nlat,file%nlev,plon,plat,pcols,ncols,file%weight_x,file%weight_y,wrk3d_in, &
- loc_arr(:,:,c-begchunk+1), lons,lats,file%count_x,file%count_y,file%index_x,file%index_y)
- enddo
+ call xy_interp(file%nlon,file%nlat,file%nlev,plon,plat,pcols,ncols, &
+ file%weight0_x,file%weight0_y,wrk3d_in,loc_arr(:,:,c-begchunk+1), &
+ lons,lats,file%count0_x,file%count0_y,file%index0_x,file%index0_y)
+ enddo
+ else
+ do c = begchunk,endchunk
+ ncols = get_ncols_p(c)
+ lons(:ncols) = lon_global_grid_ndx(:ncols,c)
+ lats(:ncols) = lat_global_grid_ndx(:ncols,c)
+ call xy_interp(file%nlon,file%nlat,file%nlev,plon,plat,pcols,ncols,&
+ file%weight_x,file%weight_y,wrk3d_in,loc_arr(:,:,c-begchunk+1), &
+ lons,lats,file%count_x,file%count_y,file%index_x,file%index_y)
+ enddo
+ endif
call t_stopf('xy_interp')
else
@@ -1839,7 +1965,7 @@ subroutine interpolate_trcdata( state, flds, file, pbuf2d )
else if(file%conserve_column) then
call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, &
datain, data_out(:,:), &
- file%p0,ps,file%hyai,file%hybi)
+ file%p0,ps,file%hyai,file%hybi,file%dist)
else
call vert_interp(ncol, file%nlev, pin, state(c)%pmid, datain, data_out(:,:) )
endif
@@ -2263,20 +2389,20 @@ subroutine interpz_conserve( nsrc, ntrg, src_x, trg_x, src, trg)
do i = 1, ntrg
tl = trg_x(i)
- if ( (tl.lt.src_x(nsrc+1)).and.(trg_x(i+1).gt.src_x(1)) ) then
+ if ( (tlsrc_x(1)) ) then
do sil = 1,nsrc
- if ( (tl-src_x(sil))*(tl-src_x(sil+1)).le.0.0_r8 ) then
+ if ( (tl-src_x(sil))*(tl-src_x(sil+1))<=0.0_r8 ) then
exit
end if
end do
- if ( tl.lt.src_x(1) ) sil = 1
+ if ( tlsrc_x(j+1) ) then
y = y+(src_x(j+1)-bot)*src(j)/(src_x(j+1)-src_x(j))
bot = src_x(j+1)
else
@@ -2290,12 +2416,12 @@ subroutine interpz_conserve( nsrc, ntrg, src_x, trg_x, src, trg)
end if
end do
- if ( trg_x(1).gt.src_x(1) ) then
+ if ( trg_x(1)>src_x(1) ) then
top = trg_x(1)
bot = src_x(1)
y = 0.0_r8
do j = 1, nsrc
- if ( top.gt.src_x(j+1) ) then
+ if ( top>src_x(j+1) ) then
y = y+(src_x(j+1)-bot)*src(j)/(src_x(j+1)-src_x(j))
bot = src_x(j+1)
else
@@ -2310,7 +2436,7 @@ subroutine interpz_conserve( nsrc, ntrg, src_x, trg_x, src, trg)
end subroutine interpz_conserve
!------------------------------------------------------------------------------
- subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, hybi)
+ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, hybi, use_flight_distance)
implicit none
@@ -2320,6 +2446,7 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai,
real(r8) :: src_x(nsrc+1) ! source coordinates
real(r8), intent(in) :: trg_x(pcols,ntrg+1) ! target coordinates
real(r8), intent(in) :: src(pcols,nsrc) ! source array
+ logical, intent(in) :: use_flight_distance ! .true. = flight distance, .false. = mixing ratio
real(r8), intent(out) :: trg(pcols,ntrg) ! target array
real(r8) :: ps(pcols), p0, hyai(nsrc+1), hybi(nsrc+1)
@@ -2335,59 +2462,77 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai,
do n = 1,ncol
- do i=1,nsrc+1
- src_x(i) = p0*hyai(i)+ps(n)*hybi(i)
- enddo
-
- do i = 1, ntrg
- tl = trg_x(n,i+1)
- if( (tl.gt.src_x(1)).and.(trg_x(n,i).lt.src_x(nsrc+1)) ) then
- do sil = 1,nsrc
- if( (tl-src_x(sil))*(tl-src_x(sil+1)).le.0.0_r8 ) then
- exit
- end if
- end do
+ do i=1,nsrc+1
+ src_x(i) = p0*hyai(i)+ps(n)*hybi(i)
+ enddo
- if( tl.gt.src_x(nsrc+1)) sil = nsrc
+ do i = 1, ntrg
+ tl = trg_x(n,i+1)
+ if( (tl>src_x(1)).and.(trg_x(n,i)src_x(nsrc+1)) sil = nsrc
+
+ y = 0.0_r8
+ bot = min(tl,src_x(nsrc+1))
+ top = trg_x(n,i)
+ do j = sil,1,-1
+ if( top pin(i,levsiz)) then
dataout(i,k) = datain(i,levsiz)
else
dpu = pmid(i,k) - pin(i,kupper(i))
diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90
index e8da28d94c..4412ae963a 100644
--- a/src/physics/cam/physpkg.F90
+++ b/src/physics/cam/physpkg.F90
@@ -2014,6 +2014,7 @@ subroutine tphysbc (ztodt, state, &
use micro_mg_cam, only: massless_droplet_destroyer
use cam_snapshot, only: cam_snapshot_all_outfld_tphysbc
use cam_snapshot, only: cam_snapshot_ptend_outfld
+ use ssatcontrail, only: ssatcontrail_d0
! Arguments
@@ -2420,6 +2421,13 @@ subroutine tphysbc (ztodt, state, &
prec_pcw_macmic = 0._r8
snow_pcw_macmic = 0._r8
+ ! contrail parameterization
+ ! see Chen et al., 2012: Global contrail coverage simulated
+ ! by CAM5 with the inventory of 2006 global aircraft emissions, JAMES
+ ! https://doi.org/10.1029/2011MS000105
+ call ssatcontrail_d0(state, pbuf, ztodt, ptend)
+ call physics_update(state, ptend, ztodt, tend)
+
! initialize ptend structures where macro and microphysics tendencies are
! accumulated over macmic substeps
call physics_ptend_init(ptend_macp_all,state%psetcols,'macrophysics',lu=.true.,lv=.true.)
diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90
new file mode 100644
index 0000000000..17c620ae24
--- /dev/null
+++ b/src/physics/cam/ssatcontrail.F90
@@ -0,0 +1,232 @@
+module ssatcontrail
+! contrail parameterization
+! see Chen et al., 2012: Global contrail coverage simulated
+! by CAM5 with the inventory of 2006 global aircraft emissions, JAMES
+! https://doi.org/10.1029/2011MS000105
+ use shr_kind_mod, only: r8 => shr_kind_r8
+ use ppgrid, only: pcols, pver
+ use physics_types, only: physics_state, physics_ptend, physics_ptend_init
+ use physconst, only: cpair,mwdry,mwh2o, gravit, zvir, rair, pi, rearth, tmelt
+ use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_old_tim_idx
+ use constituents, only: cnst_get_ind, pcnst
+ use phys_grid, only: get_wght_all_p
+ use wv_saturation, only: qsat_water, qsat_ice
+ use aircraft_emit, only: get_aircraft
+
+ implicit none
+ private
+ save
+
+ public ssatcontrail_d0
+
+ ! Private data
+ real(r8), parameter :: rhoi = 500.0_r8 ! density of ice (500 kg/m3)
+ real(r8), parameter :: radius = 3.75e-6_r8 ! diameter of ice particle = 7.5 microns
+
+
+contains
+
+ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc)
+ implicit none
+
+ type(physics_state), intent(in) :: state1
+ type(physics_ptend), intent(inout) :: ptend_loc
+ type(physics_buffer_desc), pointer :: pbuf(:)
+ real(r8), intent(in) :: dtime ! time step
+!------------------------Local storage------------------------------------------------------
+ real(r8) :: Ma, Mh2o, epsi, Q, eta, p, G, T_contr, eslTc, eslT, RH_contr, qslTc, qslT
+ real(r8) :: w, esiT, qsiT, ws, RH, ei
+ integer :: i,k
+ integer :: lchnk,ncol
+ real(r8) :: contrail(pcols,pver), pcontrail(pcols,pver)
+ real(r8), pointer, dimension(:,:) :: cld ! cloud fraction
+ real(r8), pointer, dimension(:,:) :: ac_H2O
+ real(r8), pointer, dimension(:,:) :: ac_SLANT_DIST
+ integer :: itim, ifld
+ integer :: ixcldice, ixcldliq ! indices for CLDICE and CLDLIQ
+ integer :: ixnumice, ixnumliq
+ real(r8):: zi, zm, rog
+ logical :: has_aircraft_H2O
+ logical :: has_aircraft_distance
+ real(r8) :: hkl, hkk, tv
+ real(r8) :: particle_mass
+ real(r8) :: ICIWC0(pcols,pver), ICIWC, rho
+ real(r8) :: qs
+ real(r8) :: wght(pcols)
+ real(r8) :: dz, ratio(pcols,pver)
+ real(r8) :: dcld(pcols,pver)
+ real(r8) :: ac_Q, ac_Q1, ac_Q2
+ real(r8) :: RHcts(pcols,pver)
+ logical :: lq(pcnst)
+
+ integer :: yr, mon, day, ncsec
+ real(r8) :: curr_factor
+ integer :: aircraft_cnt
+ character(len=16) :: spc_name_list(30)
+
+ has_aircraft_H2O = .false.
+ has_aircraft_distance = .false.
+
+ ! Update constituents, all schemes use time split q: no tendency kept
+ call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
+ call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
+ ! Check for number concentration of cloud liquid and cloud ice (if not present)
+ ! the indices will be set to -1)
+ call cnst_get_ind('NUMICE', ixnumice, abort=.false.)
+ call cnst_get_ind('NUMLIQ', ixnumliq, abort=.false.)
+
+ call get_aircraft(aircraft_cnt, spc_name_list)
+!-----------------------------------------------------------------------------------------
+! check if ac_H2O in namelist
+! if not, then bypass this subroutine
+!-----------------------------------------------------------------------------------------
+ if(aircraft_cnt>0) then
+ do i=1,aircraft_cnt
+ if(trim(spc_name_list(i)) == 'ac_H2O') then
+ has_aircraft_H2O = .true.
+ endif
+ if(trim(spc_name_list(i)) == 'ac_SLANT_DIST' .or. trim(spc_name_list(i)) == 'ac_TRACK_DIST') then
+ has_aircraft_distance = .true.
+ endif
+ enddo
+ endif
+!------------------------------------------------------------------------------------------
+ lq(:) = .FALSE.
+ lq(1) = .TRUE.
+ lq(ixcldice) = .TRUE.
+ lq(ixnumice) = .TRUE.
+
+ call physics_ptend_init(ptend_loc, state1%psetcols,'ssatcontrail',ls=.true.,lq = lq)
+!-----------------------------------------------------------------------------------------
+ if(.not. has_aircraft_H2O) return
+ if(.not. has_aircraft_distance) return
+
+ particle_mass = 4._r8/3._r8*pi*rhoi*radius**3 ! mass of ice particle
+
+ rog = rair/gravit ! Rd/Cp
+
+
+ contrail(:,:) = 0.0_r8
+ pcontrail(:,:) = 0.0_r8
+ ICIWC0(:,:) = 0.0_r8
+ RHcts(:,:) = 0.0_r8
+
+ lchnk = state1%lchnk
+ ncol = state1%ncol
+
+ call get_wght_all_p(lchnk, ncol, wght)
+
+ itim = pbuf_old_tim_idx()
+ ifld = pbuf_get_index('CLD')
+ call pbuf_get_field(pbuf, ifld, cld, (/1,1,itim/),(/pcols,pver,1/))
+
+ ifld = pbuf_get_index('ac_H2O')
+ call pbuf_get_field(pbuf,ifld,ac_H2O)
+ ifld = pbuf_get_index('ac_SLANT_DIST')
+ call pbuf_get_field(pbuf,ifld,ac_SLANT_DIST)
+
+! adjust h2o to volume mixing ratio (mass adjustment and conversion from g/kg to kg/kg)
+ Ma = mwdry
+ Mh2o = mwh2o
+
+! contrail paramter (G = CF*p/epi)
+! and Schumann 1996 DOI: 10.1127/metz/5/1996/4, reprinted by Ponater 2002, JGR (eq 6-8) DOI: 10.1029/2011MS000105
+
+ epsi = Mh2o/Ma
+ ei = 1.21_r8 ! water vapor emision index (g) h2o per kg fuel (Schumann 96)
+ Q = 43.e6_r8 ! specific combustion heat Schummann 1996, Q = 43 MJ/kg
+ eta = 0.3_r8 ! propulsion effieciency (Ponater 2002)
+
+ ratio = 0._r8
+ dcld = 0._r8
+
+ do i=1,ncol
+ do k=1,pver
+ p = (state1%pint(i,k)+state1%pint(i,k+1))/2.0_r8
+
+ G = (ei*cpair*p)/(epsi*Q*(1.0_r8-eta)) ! eq 7, Ponater JGR 2002
+
+ T_contr = -46.46_r8+9.43_r8*log(G-0.053_r8)+0.72_r8*log(G-0.053_r8)*log(G-0.053_r8) ! eq 8, Ponater JGR 2002
+ T_contr = T_contr + tmelt ! convert to Kelvin
+
+ ! compute saturation pressure
+ call qsat_water(T_contr, p, eslTc, qslTc)
+ call qsat_water(state1%t(i,k), p, eslT, qslT)
+
+ RH_contr = (G*(state1%t(i,k)-T_contr)+eslTc)/eslT
+ ! RH_contr ranges between 0 and 1
+ if(RH_contr>1.0_r8) RH_contr = 1.0_r8
+ if(RH_contr<0.0_r8) RH_contr = 0.0_r8
+
+ w = state1%q(i,k,1)/(1.0_r8-state1%q(i,k,1)) ! mixing ratio from specific humidity
+ call qsat_ice(state1%t(i,k), p, esiT, qsiT)
+ ws = epsi*esiT/(p-esiT) ! saturation mixing ration with respect to ice
+ qs = ws/(1.0_r8+ws)
+
+ RH = w/ws ! relative humidity with respect to ice
+ if( RH>=1.0_r8 ) RHcts(i,k) = 1.0_r8
+
+! Schumann, U. “Contrail Cirrus.” In Cirrus, edited by D. K. Lynch and others, 231–55. Oxford University Press, 2002
+! IWC(g/m3) = exp(6.97+0.103*T(C))*1e-3
+! IWC(kg/m3) = exp(6.97+0.103*T(C))*1e-6
+
+ ICIWC0(i,k) = exp(6.97_r8+0.103_r8*(state1%t(i,k)-tmelt)) ! in mg/m3
+ rho = p/(rair*state1%t(i,k))
+ ICIWC = ICIWC0(i,k)/rho*1.0e-6_r8
+
+
+! persistent contrail condition
+ if( (state1%t(i,k)RH_contr).and.(RH>1.0_r8).and.(ac_H2O(i,k)>0.0_r8) ) then
+
+! if persistent contrail, H2O emitted from aircraft turns into cloud ice
+ dz = state1%zi(i,k)-state1%zi(i,k+1)
+ ratio(i,k) = (ac_SLANT_DIST(i,k)*dtime*1.e4_r8)/(dz*rearth*rearth*wght(i))
+
+ ac_Q = min(ac_H2O(i,k)*dtime + (state1%q(i,k,1)-qs)*ratio(i,k),ratio(i,k)*ICIWC)
+ ptend_loc%q(i,k,ixcldice) = ac_Q/dtime
+
+! take out water vapor from q
+ ptend_loc%q(i,k,1) = -(ac_Q-ac_H2O(i,k)*dtime)/dtime
+
+! modify cloud fraction
+! by a prescribed ICIWC, we may deduce the new cloud fraction
+
+ cld(i,k) = min(1._r8, cld(i,k)+ac_Q/ICIWC)
+
+! modify cloud ice number concentration,
+! by assuming the particle size, the number of ice particles may be obtained
+ ptend_loc%q(i,k,ixnumice) = ac_Q/particle_mass/dtime
+
+ else
+! if not persistent contrail, just add ac_H2O to state1%q(1) (vapor phase)
+
+ ptend_loc%q(i,k,1) = ac_H2O(i,k)
+
+ endif
+
+
+ enddo
+
+! modify dry static energy if water field is added to any grid cell
+! this bypasses geopotential_t which assumes dry static energy conservation
+! water vapor added to the system is assumed to increase dry static energy
+! conservation of dry static energy by geopotential_t will lower temperature to compensate
+
+ zi = 0.0_r8
+ do k=pver,1,-1
+ hkl = state1%lnpint(i,k+1)-state1%lnpint(i,k)
+ hkk = 1._r8 - state1%pint(i,k) * hkl * state1%rpdel(i,k)
+
+ tv = state1%t(i,k) * (1._r8 + zvir*(state1%q(i,k,1)+ptend_loc%q(i,k,1)*dtime))
+
+ zm = zi + rog * tv * hkk
+ zi = zi + rog * tv * hkl
+
+ ptend_loc%s(i,k) = (cpair*state1%t(i,k)+gravit*zm + state1%phis(i) - state1%s(i,k) )/dtime
+ enddo
+
+ enddo
+
+ end subroutine ssatcontrail_d0
+
+end module ssatcontrail