Skip to content

cam6_3_014: contrail parameterization added #277

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cacraigucar merged 28 commits into from
Mar 23, 2021
Merged

cam6_3_014: contrail parameterization added #277

Show file tree
Hide file tree
Changes from 22 commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
2a77158
contrail parameterization added
chihchen24 Nov 17, 2020
65fb366
remove public of aircraft_cnt, spc_name_list
chihchen24 Jan 12, 2021
9d1a505
revised contrail code
chihchen24 Jan 12, 2021
7f88f71
revised contrail code 1/14/2021
chihchen24 Jan 14, 2021
78065ab
revision 1/15/2021
chihchen24 Jan 15, 2021
b27dc54
revision 1/15/2021 v2
chihchen24 Jan 15, 2021
3843460
revision 1/18/2021
chihchen24 Jan 18, 2021
381766f
revised 1/25/2021
chihchen24 Jan 25, 2021
43c2be5
revised 2/11/2021
chihchen24 Feb 11, 2021
088b253
Merge up to cam6_3_010
cacraigucar Mar 1, 2021
f61f53a
Add error handling and update MPI
cacraigucar Mar 2, 2021
fa7b53d
Fix errors from regression testing
cacraigucar Mar 3, 2021
70e0d73
revised 3/5/2021
chihchen24 Mar 5, 2021
04335ea
Additional changes based on reviewer comments
cacraigucar Mar 8, 2021
7202872
Add DOI
cacraigucar Mar 8, 2021
4d02c1c
Add book citation
cacraigucar Mar 8, 2021
ae9186e
revised 3/9/2021
chihchen24 Mar 9, 2021
822e263
revised 3/9/2021 v2
chihchen24 Mar 9, 2021
5ee8294
revised 3/9/2021 v3
chihchen24 Mar 9, 2021
ee72016
revised 3/9/2021 v4
chihchen24 Mar 9, 2021
1467554
Merge to cam6_3_012
cacraigucar Mar 9, 2021
49d07ae
replace problematic tab
cacraigucar Mar 9, 2021
f8aa85a
Fix some bugs and add cam regression test
cacraigucar Mar 17, 2021
3cc0e7f
Remove extra test
cacraigucar Mar 19, 2021
99baebe
Merge tag 'cam6_3_013' into cam_development_contrail
cacraigucar Mar 19, 2021
9954cd9
Merge tag 'cam6_3_013' into cam_development_contrail
cacraigucar Mar 22, 2021
b1b1a31
Fix MPI bug
cacraigucar Mar 22, 2021
3340aba
Update ChangeLog
cacraigucar Mar 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 33 additions & 7 deletions src/chemistry/utils/aircraft_emit.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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()

!------------------------------------------------------------------
Expand All @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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')
Expand All @@ -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
Expand Down
71 changes: 54 additions & 17 deletions src/chemistry/utils/horizontal_interpolate.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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
! |--------------------------------|
Expand All @@ -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


Expand All @@ -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
Expand Down
Loading