From 2a771584c69090603074c6ce5799e091ce0c24f0 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Mon, 16 Nov 2020 19:28:17 -0700 Subject: [PATCH 01/24] contrail parameterization added --- src/chemistry/utils/aircraft_emit.F90 | 24 +- .../utils/horizontal_interpolate.F90 | 211 ++++++++++- src/chemistry/utils/tracer_data.F90 | 155 +++++++- src/physics/cam/physpkg.F90 | 4 + src/physics/cam/ssatcontrail.F90 | 342 ++++++++++++++++++ 5 files changed, 721 insertions(+), 15 deletions(-) create mode 100644 src/physics/cam/ssatcontrail.F90 diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index 1403a4e264..112861237f 100644 --- a/src/chemistry/utils/aircraft_emit.F90 +++ b/src/chemistry/utils/aircraft_emit.F90 @@ -40,16 +40,16 @@ 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 ',& + 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_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'/) + .false., .false., .false., .false.,.false.,.false.,.false.,.false./) + character(len=3) :: mixtype(N_AERO) = (/'wet','wet','wet','wet','wet','wet','wet','wet','wet','wet','wet','wet','wet'/) real(r8) :: cptmp = 666.0_r8 real(r8) :: qmin = 0.0_r8 @@ -64,9 +64,10 @@ module aircraft_emit integer :: number_flds - integer :: aircraft_cnt = 0 - character(len=16) :: spc_name_list(N_AERO) + integer, public :: aircraft_cnt = 0 + character(len=16), public :: spc_name_list(N_AERO) character(len=256) :: spc_flist(N_AERO),spc_fname(N_AERO) + integer :: dist(N_AERO) contains @@ -90,6 +91,8 @@ subroutine aircraft_emit_register() !------------------------------------------------------------------ ! Return if air_specifier is blank (no aircraft data to process) !------------------------------------------------------------------ + dist(:) = 0 + aircraft_cnt = 0 if (air_specifier(1) == "") return ! count aircraft emission species used in the simulation @@ -108,6 +111,8 @@ subroutine aircraft_emit_register() call endrun('aircraft_emit_register: '//trim(spc_name)//' is not in the aircraft emission dataset') endif + if( mm<=2 ) dist(n) = 1 + aircraft_cnt = aircraft_cnt + 1 call pbuf_add_field(aero_names(mm),'physpkg',dtype_r8,(/pcols,pver/),idx) @@ -189,6 +194,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 + if( dist(m) == 1 ) forcings_air(m)%file%dist = .true. 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 +320,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 +347,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..7814da8c48 100644 --- a/src/chemistry/utils/horizontal_interpolate.F90 +++ b/src/chemistry/utils/horizontal_interpolate.F90 @@ -14,7 +14,7 @@ module horizontal_interpolate real(r8) :: gw1(1000), gw2(1000) - public :: xy_interp_init, xy_interp + public :: xy_interp_init, xy_interp0_init, xy_interp contains subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) @@ -202,7 +202,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) ! |---------------------------------| ! 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 + elseif ( (y1_north.gt.y2_north).and.(y1_south.lt.y2_south) ) then ! case 4: ! y1_south y1_north ! |--------------------------------| @@ -215,6 +215,213 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) enddo end subroutine xy_interp_init + subroutine xy_interp0_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) +!------------------------------------------------------------------------------------------------------------ +! 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 +! weight_y(jm2,jm1) is the weighting function for meridional interpolation +! +! Author: Chih-Chieh (Jack) Chen -- May 2010 +! +!------------------------------------------------------------------------------------------------------------ + implicit none + integer, intent(in) :: im1, jm1, im2, jm2 + real(r8), intent(in) :: lon0(im1), lat0(jm1) + real(r8), intent(out) :: weight_x(im2,im1), weight_y(jm2,jm1) + + real(r8) :: lon1(im1), lat1(jm1) + real(r8) :: lon2(im2), lat2(jm2) + real(r8) :: slon1(im1+1), slon2(im2+1), slat1(jm1+1), slat2(jm2+1) + real(r8) :: x1_west, x1_east, x2_west, x2_east + real(r8) :: y1_south, y1_north, y2_south, y2_north + integer :: i1, j1, i2, j2 + + weight_x(:,:) = 0.0_r8 + weight_y(:,:) = 0.0_r8 + +! lon0 & lat0 are longitude & latitude on the source mesh in radians +! convert lon1, lat1 from radians to degrees + lon1(:) = lon0(:)/SHR_CONST_PI*180.0_r8 + lat1(:) = lat0(:)/SHR_CONST_PI*180.0_r8 + +! set up lon2, lat2 (target mesh), in CAM convention + do i2=1,im2 + lon2(i2) = (float(i2)-1.0_r8)*360.0_r8/float(im2) + enddo + do j2=1,jm2 + lat2(j2) = -90.0_r8+(float(j2)-1.0_r8)*180.0_r8/(float(jm2)-1.0_r8) + enddo + + +! set up staggered longitudes (cell edges in x) + do i1=2,im1 + slon1(i1) = (lon1(i1-1)+lon1(i1))/2.0_r8 + enddo + slon1(1) = lon1(1)-(lon1(2)-lon1(1))/2.0_r8 + slon1(im1+1) = lon1(im1)+(lon1(im1)-lon1(im1-1))/2.0_r8 + + do i2=2,im2 + slon2(i2) = (lon2(i2-1)+lon2(i2))/2.0_r8 + enddo + slon2(1) = lon2(1)-(lon2(2)-lon2(1))/2.0_r8 + slon2(im2+1) = lon2(im2)+(lon2(im2)-lon2(im2-1))/2.0_r8 + +! set up staggered lattiudes (cell edges in y) + slat1(1)=-90.0_r8 + do j1=2,jm1 + slat1(j1) = (lat1(j1-1)+lat1(j1))/2.0_r8 + enddo + slat1(jm1+1)=90.0_r8 + + slat2(1)=-90.0_r8 + do j2=2,jm2 + slat2(j2)=(lat2(j2-1)+lat2(j2))/2.0_r8 + enddo + slat2(jm2+1)=90.0_r8 + +! compute Guassian weight for two meshes (discrete form of cos(lat).) + do j1=1,jm1 + gw1(j1) = sin(slat1(j1+1)/180.0_r8*SHR_CONST_PI)-sin(slat1(j1)/180.0_r8*SHR_CONST_PI) + enddo + + do j2=1,jm2 + gw2(j2) = sin(slat2(j2+1)/180.0_r8*SHR_CONST_PI)-sin(slat2(j2)/180.0_r8*SHR_CONST_PI) + enddo + + +! add 360 to slon1 and slon2 + slon1(:) = slon1(:)+360.0_r8 + slon2(:) = slon2(:)+360.0_r8 + + do i2=1,im2 + +! target grid east-west boundaries + x2_west=slon2(i2) + x2_east=slon2(i2+1) + + do i1=1,im1 + +! source grid east-west boundaries + x1_west=slon1(i1) + x1_east=slon1(i1+1) + +! 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 +! case 1: +! x1_west x1_east +! |-------------------| +! |---------------------------------| +! x2_west x2_east +! weight_x(i2,i1) = (x1_east-x1_west)/(x2_east-x2_west) + weight_x(i2,i1) = 1.0_r8 + elseif ( (x1_west.ge.x2_west).and.(x1_west.lt.x2_east) ) then +! case 2: +! x1_west x1_east +! |--------------------------------| +! |---------------------------------| +! x2_west x2_east +! weight_x(i2,i1) = (x2_east-x1_west)/(x2_east-x2_west) + weight_x(i2,i1) = (x2_east-x1_west)/(x1_east-x1_west) + elseif ( (x1_east>x2_west).and.(x1_east.le.x2_east) ) then +! case 3: +! x1_west x1_east +! |--------------------------------| +! |---------------------------------| +! x2_west x2_east +! weight_x(i2,i1) = (x1_east-x2_west)/(x2_east-x2_west) + weight_x(i2,i1) = (x1_east-x2_west)/(x1_east-x1_west) + elseif ( (x1_east.gt.x2_east).and.(x1_west.lt.x2_west) ) then +! case 4: +! x1_west x1_east +! |--------------------------------| +! |--------------------| +! x2_west x2_east + weight_x(i2,i1) = 1._r8 + endif + + enddo + enddo + + +! consider end points + if(slon1(im1+1).gt.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)) + weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon1(im1+1)-slon1(im1)) + endif + + if(slon1(im1+1).lt.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)) + weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon1(2)-slon1(1)) + endif + + + + do j2=1,jm2 +! target grid north-south boundaries + y2_south=slat2(j2) + y2_north=slat2(j2+1) + + do j1=1,jm1 + +! source grid north-south boundaries + y1_south=slat1(j1) + y1_north=slat1(j1+1) + +! 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 +! note: there is Guassian weight to consider in the meridional direction! + + if( (y1_south.ge.y2_south).and.(y1_north.le.y2_north) ) then +! case 1: +! y1_south y1_north +! |-------------------| +! |---------------------------------| +! y2_south y2_north +! weight_y(j2,j1) = gw1(j1)/gw2(j2) + weight_y(j2,j1) = 1.0_r8 + elseif ( (y1_south.ge.y2_south).and.(y1_south.lt.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) + weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south) + elseif ( (y1_north.gt.y2_south).and.(y1_north.le.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) + weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south) + elseif ( (y1_north.gt.y2_north).and.(y1_south.lt.y2_south) ) then +! case 4: +! y1_south y1_north +! |--------------------------------| +! |---------------------| +! y2_south y2_north + weight_y(j2,j1) = 1._r8 + endif + + enddo + enddo + + end subroutine xy_interp0_init + subroutine xy_interp(im1,jm1,km1,im2,jm2,pcols,ncols,weight_x,weight_y,var_src,var_trg,lons,lats,count_x,count_y,index_x,index_y) !------------------------------------------------------------------------------------------------------------- diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 64ad69dd99..3a78ed3638 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -110,6 +110,12 @@ module tracer_data real(r8), pointer, dimension(:,:) :: weight_x => 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 = .false. + real(r8) :: p0 type(var_desc_t) :: ps_id logical, allocatable, dimension(:) :: in_pbuf @@ -161,7 +167,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & use dyn_grid, only : get_dyn_grid_parm, get_horiz_grid_d 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 + use horizontal_interpolate, only : xy_interp_init, xy_interp0_init #if ( defined SPMD ) use mpishorthand, only: mpicom, mpir8, mpiint #endif @@ -588,8 +594,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) @@ -630,6 +636,19 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & file%index_x(:,:) = 0 file%index_y(:,:) = 0 + allocate(file%weight0_x(plon,file%nlon)) + allocate(file%weight0_y(plat,file%nlat)) + allocate(file%count0_x(plon)) + allocate(file%count0_y(plat)) + allocate(file%index0_x(plon,file%nlon)) + allocate(file%index0_y(plat,file%nlat)) + 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 + 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) @@ -653,6 +672,28 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & endif enddo enddo + + call xy_interp0_init(file%nlon,file%nlat,file%lons,file%lats,plon,plat,file%weight0_x,file%weight0_y) + + do i2=1,plon + file%count0_x(i2) = 0 + do i1=1,file%nlon + if(file%weight0_x(i2,i1).gt.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).gt.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 #if ( defined SPMD) @@ -662,6 +703,12 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & 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) + call mpibcast(file%weight0_x, plon*file%nlon, mpir8 , 0, mpicom) + call mpibcast(file%weight0_y, plat*file%nlat, mpir8 , 0, mpicom) + call mpibcast(file%count0_x, plon, mpiint , 0, mpicom) + call mpibcast(file%count0_y, plat, mpiint , 0, mpicom) + call mpibcast(file%index0_x, plon*file%nlon, mpiint , 0, mpicom) + call mpibcast(file%index0_y, plat*file%nlat, mpiint , 0, mpicom) #endif endif @@ -1623,16 +1670,25 @@ subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order) if(file%weight_by_lat) then call t_startf('xy_interp') + 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) - do c = begchunk,endchunk + 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 - + enddo + endif call t_stopf('xy_interp') else @@ -1834,9 +1890,15 @@ subroutine interpolate_trcdata( state, flds, file, pbuf2d ) if ( file%top_bndry ) then call vert_interp_ub(ncol, file%nlev, file%levs, datain(:ncol,:), data_out(:ncol,:) ) else if(file%conserve_column) then + if( file%dist ) then + call vert_rebin(ncol,file%nlev,pver,state(c)%pint, & + datain,data_out(:,:), & + file%p0,ps,file%hyai,file%hybi) + else call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, & datain, data_out(:,:), & file%p0,ps,file%hyai,file%hybi) + endif else call vert_interp(ncol, file%nlev, pin, state(c)%pmid, datain, data_out(:,:) ) endif @@ -2389,6 +2451,87 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, enddo end subroutine vert_interp_mixrat +!------------------------------------------------------------------------------ + subroutine vert_rebin( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, hybi) + + implicit none + + integer, intent(in) :: ncol + integer, intent(in) :: nsrc ! dimension source array + integer, intent(in) :: ntrg ! dimension target array + 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 + real(r8), intent(out) :: trg(pcols,ntrg) ! target array + + real(r8) :: ps(pcols), p0, hyai(nsrc+1), hybi(nsrc+1) + !--------------------------------------------------------------- + ! ... local variables + !--------------------------------------------------------------- + integer :: i, j, n + integer :: sil + real(r8) :: tl, y + real(r8) :: bot, top + + 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 + + if( tl.gt.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.lt.src_x(j) ) then + y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) + bot = src_x(j) + else + y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) + exit + endif + enddo + trg(n,i) = y + else + trg(n,i) = 0.0_r8 + end if + end do + + if( trg_x(n,ntrg+1).lt.src_x(nsrc+1) ) then + top = trg_x(n,ntrg+1) + bot = src_x(nsrc+1) + y = 0.0_r8 + do j=nsrc,1,-1 + if( top.lt.src_x(j) ) then + y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) + bot = src_x(j) + else + y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) + exit + endif + enddo + trg(n,ntrg) = trg(n,ntrg)+y + endif + +! turn mass into mixing ratio +! do i=1,ntrg +! trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i)) +! enddo + + enddo + + end subroutine vert_rebin !------------------------------------------------------------------------------ subroutine vert_interp( ncol, levsiz, pin, pmid, datain, dataout ) !-------------------------------------------------------------------------- diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 38cc4fca9d..2370cd416d 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -1889,6 +1889,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 @@ -2273,6 +2274,9 @@ subroutine tphysbc (ztodt, state, & prec_pcw_macmic = 0._r8 snow_pcw_macmic = 0._r8 + call ssatcontrail_d0(state, pbuf, ztodt, ptend, tend) + call physics_update(state, ptend, ztodt, tend) + do macmic_it = 1, cld_macmic_num_steps !=================================================== diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 new file mode 100644 index 0000000000..2ca7096df6 --- /dev/null +++ b/src/physics/cam/ssatcontrail.F90 @@ -0,0 +1,342 @@ +module ssatcontrail + + use shr_kind_mod, only: r8 => shr_kind_r8 + use ppgrid, only: pcols, pver + use cam_history, only: outfld + use cam_logfile, only: iulog + use physics_types, only: physics_state, physics_ptend, physics_tend + use physics_types, only: physics_ptend_sum, physics_update + use physics_types, only: physics_state_copy, physics_ptend_init + use physconst, only: cpair,mwdry,mwh2o, gravit, zvir, rair, pi, rearth +! use phys_buffer, only: pbuf_size_max, pbuf_fld, pbuf_old_tim_idx, pbuf_get_fld_idx, pbuf_times + use physics_buffer, only : pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_set_field, pbuf_old_tim_idx + use constituents, only: cnst_get_ind, pcnst + use aircraft_emit, only: aircraft_cnt, spc_name_list + use geopotential, only: geopotential_dse + use phys_grid, only : get_wght_all_p + use time_manager, only: get_curr_date + + implicit none + private + save + + public ssatcontrail_d0 + + ! Private data +! real(r8), parameter :: ICIWC0 = 2.0e-6 ! ICIWC = 2 mg/m3, converted to kg/kg here by dividing rho_air + real(r8), parameter :: rhoi = 500.0_r8 ! density of ice (500 kg/m3) + real(r8), parameter :: radius = 3.75e-6 ! diameter of ice particle = 7.5 microns + + +contains + + subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) + implicit none + + type(physics_state), intent(inout) :: state1 + type(physics_ptend), intent(inout) :: ptend_loc + type(physics_tend) :: tend + 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 + real(r8) :: w, esiT, 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 + 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) + integer :: itype + logical :: lq(pcnst) + + integer :: yr, mon, day, ncsec + real(r8) :: curr_factor + +! ICIWC = ICIWC0/rhodair + + has_aircraft_H2O = .false. + +!----------------------------------------------------------------------------------------- +! 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 + enddo + endif + if(.not. has_aircraft_H2O) 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) + + + ! 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.) + +!------------------------------------------------------------------------------------------ + lq(:) = .FALSE. + lq(1) = .TRUE. + lq(ixcldice) = .TRUE. + lq(ixnumice) = .TRUE. + + call physics_ptend_init(ptend_loc, state1%psetcols,'ssatcontrail',ls=.true.,lq = lq) +!----------------------------------------------------------------------------------------- + +! 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, reprinted by Ponater 2002, JGR (eq 6-8) + + epsi = Mh2o/Ma + ei = 1.21_r8 ! water vapor emiision index (g) h2o per kg fuel (Schumann 96)? + Q = 43.e6 ! specific combustion heat Schummann 1996, Q = 43 MJ/kg + eta = 0.3_r8 ! propulsion effieciency (Ponater 2002) + + ratio(i,k) = 0._r8 + dcld(i,k) = 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 + 273.15_r8 ! convert to Kelvin + + ! compute saturation pressure + itype = 0 + call gffgch(T_contr,eslTc,itype) + itype = 0 + call gffgch(state1%t(i,k),eslT,itype) + + RH_contr = (G*(state1%t(i,k)-T_contr)+eslTc)/eslT + ! RH_contr ranges between 0 and 1 + if(RH_contr.gt.1.0_r8) RH_contr = 1.0_r8 + if(RH_contr.lt.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 + itype = 1 + call gffgch(state1%t(i,k),esiT,itype) ! saturation water vapor with respect to ice + 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.ge.1.0_r8 ) RHcts(i,k) = 1.0_r8 + +! Schumann, 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)-273.16_r8)) ! in mg/m3 + rho = p/(287._r8*state1%t(i,k)) + ICIWC = ICIWC0(i,k)/rho*1.0e-6 + + +! persistent contrail condition + if( (state1%t(i,k).lt.T_contr).and.(RH.gt.RH_contr).and.(RH.gt.1.0_r8).and.(ac_H2O(i,k).gt.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 + + call physics_update(state1, ptend_loc, dtime, tend) + + end subroutine ssatcontrail_d0 + +subroutine gffgch(t ,es ,itype ) + use shr_kind_mod, only: r8 => shr_kind_r8 + use physconst, only: tmelt + use abortutils, only: endrun + use cam_logfile, only: iulog + + implicit none +!------------------------------Arguments-------------------------------- +! +! Input arguments +! + real(r8), intent(in) :: t ! Temperature +! +! Output arguments +! + integer, intent(inout) :: itype ! Flag for ice phase and associated transition + + real(r8), intent(out) :: es ! Saturation vapor pressure +! +!---------------------------Local variables----------------------------- + real(r8) e1 ! Intermediate scratch variable for es over water + real(r8) e2 ! Intermediate scratch variable for es over water + real(r8) eswtr ! Saturation vapor pressure over water + real(r8) f ! Intermediate scratch variable for es over water + real(r8) f1 ! Intermediate scratch variable for es over water + real(r8) f2 ! Intermediate scratch variable for es over water + real(r8) f3 ! Intermediate scratch variable for es over water + real(r8) f4 ! Intermediate scratch variable for es over water + real(r8) f5 ! Intermediate scratch variable for es over water + real(r8) ps ! Reference pressure (mb) + real(r8) t0 ! Reference temperature (freezing point of water) + real(r8) term1 ! Intermediate scratch variable for es over ice + real(r8) term2 ! Intermediate scratch variable for es over ice + real(r8) term3 ! Intermediate scratch variable for es over ice + real(r8) tr ! Transition range for es over water to es over ice + real(r8) ts ! Reference temperature (boiling point of water) + real(r8) weight ! Intermediate scratch variable for es transition + integer itypo ! Intermediate scratch variable for holding itype +! +!----------------------------------------------------------------------- +! +! Check on whether there is to be a transition region for es +! + if (itype < 0) then + tr = abs(real(itype,r8)) + itypo = itype + itype = 1 + else + tr = 0.0_r8 + itypo = itype + end if + if (tr > 40.0_r8) then + write(iulog,900) tr + call endrun ('GFFGCH') ! Abnormal termination + end if +! + if(t < (tmelt - tr) .and. itype == 1) go to 10 +! +! Water +! + ps = 1013.246_r8 + ts = 373.16_r8 + e1 = 11.344_r8*(1.0_r8 - t/ts) + e2 = -3.49149_r8*(ts/t - 1.0_r8) + f1 = -7.90298_r8*(ts/t - 1.0_r8) + f2 = 5.02808_r8*log10(ts/t) + f3 = -1.3816_r8*(10.0_r8**e1 - 1.0_r8)/10000000.0_r8 + f4 = 8.1328_r8*(10.0_r8**e2 - 1.0_r8)/1000.0_r8 + f5 = log10(ps) + f = f1 + f2 + f3 + f4 + f5 + es = (10.0_r8**f)*100.0_r8 + eswtr = es +! + if(t >= tmelt .or. itype == 0) go to 20 +! +! Ice +! +10 continue + t0 = tmelt + term1 = 2.01889049_r8/(t0/t) + term2 = 3.56654_r8*log(t0/t) + term3 = 20.947031_r8*(t0/t) + es = 575.185606e10_r8*exp(-(term1 + term2 + term3)) +! + if (t < (tmelt - tr)) go to 20 +! +! Weighted transition between water and ice +! + weight = min((tmelt - t)/tr,1.0_r8) + es = weight*es + (1.0_r8 - weight)*eswtr +! +20 continue + itype = itypo + return +! +900 format('GFFGCH: FATAL ERROR ******************************',/, & + 'TRANSITION RANGE FOR WATER TO ICE SATURATION VAPOR', & + ' PRESSURE, TR, EXCEEDS MAXIMUM ALLOWABLE VALUE OF', & + ' 40.0 DEGREES C',/, ' TR = ',f7.2) +! +end subroutine gffgch + + +end module ssatcontrail From 65fb3664cba65e6a09d043b5837f79e775c42a05 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Mon, 11 Jan 2021 22:18:38 -0700 Subject: [PATCH 02/24] remove public of aircraft_cnt, spc_name_list --- src/chemistry/utils/aircraft_emit.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index 112861237f..5ef577101c 100644 --- a/src/chemistry/utils/aircraft_emit.F90 +++ b/src/chemistry/utils/aircraft_emit.F90 @@ -64,8 +64,8 @@ module aircraft_emit integer :: number_flds - integer, public :: aircraft_cnt = 0 - character(len=16), public :: spc_name_list(N_AERO) + integer :: aircraft_cnt = 0 + character(len=16) :: spc_name_list(N_AERO) character(len=256) :: spc_flist(N_AERO),spc_fname(N_AERO) integer :: dist(N_AERO) From 9d1a505591544606cf16ac57833b9dab7523f7af Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Tue, 12 Jan 2021 00:36:40 -0700 Subject: [PATCH 03/24] revised contrail code --- src/chemistry/utils/aircraft_emit.F90 | 17 +- .../utils/horizontal_interpolate.F90 | 249 +++--------------- src/chemistry/utils/tracer_data.F90 | 137 +++------- src/physics/cam/physpkg.F90 | 3 + src/physics/cam/ssatcontrail.F90 | 134 ++-------- 5 files changed, 115 insertions(+), 425 deletions(-) diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index 5ef577101c..ea05fa58d8 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 @@ -71,6 +72,20 @@ module aircraft_emit 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 + + cnt = aircraft_cnt + if( cnt.gt. 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() !------------------------------------------------------------------ @@ -111,7 +126,7 @@ subroutine aircraft_emit_register() call endrun('aircraft_emit_register: '//trim(spc_name)//' is not in the aircraft emission dataset') endif - if( mm<=2 ) dist(n) = 1 + if (trim(spc_name) == 'ac_SLANT_DIST'.or. trim(spc_name) == 'ac_TRACK_DIST') dist(n) = 1 aircraft_cnt = aircraft_cnt + 1 call pbuf_add_field(aero_names(mm),'physpkg',dtype_r8,(/pcols,pver/),idx) diff --git a/src/chemistry/utils/horizontal_interpolate.F90 b/src/chemistry/utils/horizontal_interpolate.F90 index 7814da8c48..c2d3212b64 100644 --- a/src/chemistry/utils/horizontal_interpolate.F90 +++ b/src/chemistry/utils/horizontal_interpolate.F90 @@ -14,10 +14,10 @@ module horizontal_interpolate real(r8) :: gw1(1000), gw2(1000) - public :: xy_interp_init, xy_interp0_init, xy_interp + 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,ii) !------------------------------------------------------------------------------------------------------------ ! 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 @@ -27,7 +27,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) ! !------------------------------------------------------------------------------------------------------------ implicit none - integer, intent(in) :: im1, jm1, im2, jm2 + integer, intent(in) :: im1, jm1, im2, jm2, ii real(r8), intent(in) :: lon0(im1), lat0(jm1) real(r8), intent(out) :: weight_x(im2,im1), weight_y(jm2,jm1) @@ -116,21 +116,33 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) ! |-------------------| ! |---------------------------------| ! x2_west x2_east + if(ii.eq.0) then + weight_x(i2,i1) = 1.0_r8 + else weight_x(i2,i1) = (x1_east-x1_west)/(x2_east-x2_west) + endif elseif ( (x1_west.ge.x2_west).and.(x1_west.lt.x2_east) ) then ! case 2: ! x1_west x1_east ! |--------------------------------| ! |---------------------------------| ! x2_west x2_east + if(ii.eq.0) 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) + endif elseif ( (x1_east>x2_west).and.(x1_east.le.x2_east) ) then ! case 3: ! x1_west x1_east ! |--------------------------------| ! |---------------------------------| ! x2_west x2_east + if(ii.eq.0) 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) + endif elseif ( (x1_east.gt.x2_east).and.(x1_west.lt.x2_west) ) then ! case 4: ! x1_west x1_east @@ -151,7 +163,11 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) ! |-------------------------| ! |----------------|......................| ! slon2(im2) slon2(im2+1) slon2(2) (note: slon2(im2+1) = slon2(1)) + if(ii.eq.0) 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 @@ -160,7 +176,11 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) ! |-------------------------|.............................| ! |-------------------------------| ! slon2(im2) slon2(im2+1) <--- end point + if(ii.eq.0) 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 @@ -187,241 +207,50 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) ! |-------------------| ! |---------------------------------| ! y2_south y2_north + if(ii.eq.0) then + weight_y(j2,j1) = 1.0_r8 + else weight_y(j2,j1) = gw1(j1)/gw2(j2) + endif elseif ( (y1_south.ge.y2_south).and.(y1_south.lt.y2_north) ) then ! case 2: ! y1_south y1_north ! |--------------------------------| ! |---------------------------------| ! y2_south y2_north + if(ii.eq.0) 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.gt.y2_south).and.(y1_north.le.y2_north) ) then ! case 3: ! y1_south y1_north ! |--------------------------------| ! |---------------------------------| ! y2_south y2_north + if(ii.eq.0) 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.gt.y2_north).and.(y1_south.lt.y2_south) ) then ! case 4: ! y1_south y1_north ! |--------------------------------| ! |---------------------| ! y2_south y2_north + if(ii.eq.0) then + weight_y(j2,j1) = 1._r8 + else weight_y(j2,j1) = gw1(j1)/gw2(j2) + endif endif enddo enddo end subroutine xy_interp_init - subroutine xy_interp0_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y) -!------------------------------------------------------------------------------------------------------------ -! 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 -! weight_y(jm2,jm1) is the weighting function for meridional interpolation -! -! Author: Chih-Chieh (Jack) Chen -- May 2010 -! -!------------------------------------------------------------------------------------------------------------ - implicit none - integer, intent(in) :: im1, jm1, im2, jm2 - real(r8), intent(in) :: lon0(im1), lat0(jm1) - real(r8), intent(out) :: weight_x(im2,im1), weight_y(jm2,jm1) - - real(r8) :: lon1(im1), lat1(jm1) - real(r8) :: lon2(im2), lat2(jm2) - real(r8) :: slon1(im1+1), slon2(im2+1), slat1(jm1+1), slat2(jm2+1) - real(r8) :: x1_west, x1_east, x2_west, x2_east - real(r8) :: y1_south, y1_north, y2_south, y2_north - integer :: i1, j1, i2, j2 - - weight_x(:,:) = 0.0_r8 - weight_y(:,:) = 0.0_r8 - -! lon0 & lat0 are longitude & latitude on the source mesh in radians -! convert lon1, lat1 from radians to degrees - lon1(:) = lon0(:)/SHR_CONST_PI*180.0_r8 - lat1(:) = lat0(:)/SHR_CONST_PI*180.0_r8 - -! set up lon2, lat2 (target mesh), in CAM convention - do i2=1,im2 - lon2(i2) = (float(i2)-1.0_r8)*360.0_r8/float(im2) - enddo - do j2=1,jm2 - lat2(j2) = -90.0_r8+(float(j2)-1.0_r8)*180.0_r8/(float(jm2)-1.0_r8) - enddo - - -! set up staggered longitudes (cell edges in x) - do i1=2,im1 - slon1(i1) = (lon1(i1-1)+lon1(i1))/2.0_r8 - enddo - slon1(1) = lon1(1)-(lon1(2)-lon1(1))/2.0_r8 - slon1(im1+1) = lon1(im1)+(lon1(im1)-lon1(im1-1))/2.0_r8 - - do i2=2,im2 - slon2(i2) = (lon2(i2-1)+lon2(i2))/2.0_r8 - enddo - slon2(1) = lon2(1)-(lon2(2)-lon2(1))/2.0_r8 - slon2(im2+1) = lon2(im2)+(lon2(im2)-lon2(im2-1))/2.0_r8 - -! set up staggered lattiudes (cell edges in y) - slat1(1)=-90.0_r8 - do j1=2,jm1 - slat1(j1) = (lat1(j1-1)+lat1(j1))/2.0_r8 - enddo - slat1(jm1+1)=90.0_r8 - - slat2(1)=-90.0_r8 - do j2=2,jm2 - slat2(j2)=(lat2(j2-1)+lat2(j2))/2.0_r8 - enddo - slat2(jm2+1)=90.0_r8 - -! compute Guassian weight for two meshes (discrete form of cos(lat).) - do j1=1,jm1 - gw1(j1) = sin(slat1(j1+1)/180.0_r8*SHR_CONST_PI)-sin(slat1(j1)/180.0_r8*SHR_CONST_PI) - enddo - - do j2=1,jm2 - gw2(j2) = sin(slat2(j2+1)/180.0_r8*SHR_CONST_PI)-sin(slat2(j2)/180.0_r8*SHR_CONST_PI) - enddo - - -! add 360 to slon1 and slon2 - slon1(:) = slon1(:)+360.0_r8 - slon2(:) = slon2(:)+360.0_r8 - - do i2=1,im2 - -! target grid east-west boundaries - x2_west=slon2(i2) - x2_east=slon2(i2+1) - - do i1=1,im1 - -! source grid east-west boundaries - x1_west=slon1(i1) - x1_east=slon1(i1+1) - -! 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 -! case 1: -! x1_west x1_east -! |-------------------| -! |---------------------------------| -! x2_west x2_east -! weight_x(i2,i1) = (x1_east-x1_west)/(x2_east-x2_west) - weight_x(i2,i1) = 1.0_r8 - elseif ( (x1_west.ge.x2_west).and.(x1_west.lt.x2_east) ) then -! case 2: -! x1_west x1_east -! |--------------------------------| -! |---------------------------------| -! x2_west x2_east -! weight_x(i2,i1) = (x2_east-x1_west)/(x2_east-x2_west) - weight_x(i2,i1) = (x2_east-x1_west)/(x1_east-x1_west) - elseif ( (x1_east>x2_west).and.(x1_east.le.x2_east) ) then -! case 3: -! x1_west x1_east -! |--------------------------------| -! |---------------------------------| -! x2_west x2_east -! weight_x(i2,i1) = (x1_east-x2_west)/(x2_east-x2_west) - weight_x(i2,i1) = (x1_east-x2_west)/(x1_east-x1_west) - elseif ( (x1_east.gt.x2_east).and.(x1_west.lt.x2_west) ) then -! case 4: -! x1_west x1_east -! |--------------------------------| -! |--------------------| -! x2_west x2_east - weight_x(i2,i1) = 1._r8 - endif - - enddo - enddo - - -! consider end points - if(slon1(im1+1).gt.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)) - weight_x(1,im1)= weight_x(1,im1)+(slon1(im1+1)-slon2(im2+1))/(slon1(im1+1)-slon1(im1)) - endif - - if(slon1(im1+1).lt.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)) - weight_x(im2,1) = weight_x(im2,1)+(slon2(1)-slon1(1))/(slon1(2)-slon1(1)) - endif - - - - do j2=1,jm2 -! target grid north-south boundaries - y2_south=slat2(j2) - y2_north=slat2(j2+1) - - do j1=1,jm1 - -! source grid north-south boundaries - y1_south=slat1(j1) - y1_north=slat1(j1+1) - -! 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 -! note: there is Guassian weight to consider in the meridional direction! - - if( (y1_south.ge.y2_south).and.(y1_north.le.y2_north) ) then -! case 1: -! y1_south y1_north -! |-------------------| -! |---------------------------------| -! y2_south y2_north -! weight_y(j2,j1) = gw1(j1)/gw2(j2) - weight_y(j2,j1) = 1.0_r8 - elseif ( (y1_south.ge.y2_south).and.(y1_south.lt.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) - weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south) - elseif ( (y1_north.gt.y2_south).and.(y1_north.le.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) - weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south) - elseif ( (y1_north.gt.y2_north).and.(y1_south.lt.y2_south) ) then -! case 4: -! y1_south y1_north -! |--------------------------------| -! |---------------------| -! y2_south y2_north - weight_y(j2,j1) = 1._r8 - endif - - enddo - enddo - - end subroutine xy_interp0_init - subroutine xy_interp(im1,jm1,km1,im2,jm2,pcols,ncols,weight_x,weight_y,var_src,var_trg,lons,lats,count_x,count_y,index_x,index_y) !------------------------------------------------------------------------------------------------------------- diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 3a78ed3638..e72429c253 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -167,10 +167,11 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & use dyn_grid, only : get_dyn_grid_parm, get_horiz_grid_d 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, xy_interp0_init + use horizontal_interpolate, only : xy_interp_init #if ( defined SPMD ) use mpishorthand, only: mpicom, mpir8, mpiint #endif + use aircraf_emit, only: get_aircraft implicit none @@ -203,6 +204,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 :: aircrat_cnt call specify_fields( specifier, flds ) @@ -594,8 +596,11 @@ 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) - 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)) + call get_aircraft(aircraft_cnt) + if(aircraft_cnt.gt.0) then + 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)) + endif lon_global_grid_ndx=-huge(1) lat_global_grid_ndx=-huge(1) @@ -651,7 +656,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & 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,1) do i2=1,plon file%count_x(i2) = 0 @@ -673,7 +678,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & enddo enddo - call xy_interp0_init(file%nlon,file%nlat,file%lons,file%lats,plon,plat,file%weight0_x,file%weight0_y) + call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats,plon,plat,file%weight0_x,file%weight0_y,0) do i2=1,plon file%count0_x(i2) = 0 @@ -1891,13 +1896,13 @@ subroutine interpolate_trcdata( state, flds, file, pbuf2d ) call vert_interp_ub(ncol, file%nlev, file%levs, datain(:ncol,:), data_out(:ncol,:) ) else if(file%conserve_column) then if( file%dist ) then - call vert_rebin(ncol,file%nlev,pver,state(c)%pint, & - datain,data_out(:,:), & - file%p0,ps,file%hyai,file%hybi) + call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, & + datain, data_out(:,:), & + file%p0,ps,file%hyai,file%hybi,0) else 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,1) endif else call vert_interp(ncol, file%nlev, pin, state(c)%pmid, datain, data_out(:,:) ) @@ -2369,7 +2374,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, ii) implicit none @@ -2379,6 +2384,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 + integer, intent(in) :: ii ! 0: rebin only, otherwise mixing ratio real(r8), intent(out) :: trg(pcols,ntrg) ! target array real(r8) :: ps(pcols), p0, hyai(nsrc+1), hybi(nsrc+1) @@ -2414,10 +2420,18 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, top = trg_x(n,i) do j = sil,1,-1 if( top.lt.src_x(j) ) then + if(ii.eq.0) then + y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) + else y = y+(bot-src_x(j))*src(n,j) + endif bot = src_x(j) else - y = y+(bot-top)*src(n,j) + if(ii.eq.0) then + y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) + else + y = y+(bot-top)*src(n,j) + endif exit endif enddo @@ -2433,10 +2447,18 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, y = 0.0_r8 do j=nsrc,1,-1 if( top.lt.src_x(j) ) then - y = y+(bot-src_x(j))*src(n,j) + if(ii.eq.0) then + y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) + else + y = y+(bot-src_x(j))*src(n,j) + endif bot = src_x(j) else - y = y+(bot-top)*src(n,j) + if(ii.eq.0) then + y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) + else + y = y+(bot-top)*src(n,j) + endif exit endif enddo @@ -2444,94 +2466,15 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, endif ! turn mass into mixing ratio - do i=1,ntrg - trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i)) - enddo - - enddo - - end subroutine vert_interp_mixrat -!------------------------------------------------------------------------------ - subroutine vert_rebin( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, hybi) - - implicit none - - integer, intent(in) :: ncol - integer, intent(in) :: nsrc ! dimension source array - integer, intent(in) :: ntrg ! dimension target array - 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 - real(r8), intent(out) :: trg(pcols,ntrg) ! target array - - real(r8) :: ps(pcols), p0, hyai(nsrc+1), hybi(nsrc+1) - !--------------------------------------------------------------- - ! ... local variables - !--------------------------------------------------------------- - integer :: i, j, n - integer :: sil - real(r8) :: tl, y - real(r8) :: bot, top - - 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 - - if( tl.gt.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.lt.src_x(j) ) then - y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) - bot = src_x(j) - else - y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) - exit - endif - enddo - trg(n,i) = y - else - trg(n,i) = 0.0_r8 - end if - end do - - if( trg_x(n,ntrg+1).lt.src_x(nsrc+1) ) then - top = trg_x(n,ntrg+1) - bot = src_x(nsrc+1) - y = 0.0_r8 - do j=nsrc,1,-1 - if( top.lt.src_x(j) ) then - y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) - bot = src_x(j) - else - y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) - exit - endif + if(ii.ne.0) then + do i=1,ntrg + trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i)) enddo - trg(n,ntrg) = trg(n,ntrg)+y endif - -! turn mass into mixing ratio -! do i=1,ntrg -! trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i)) -! enddo - + enddo - end subroutine vert_rebin + end subroutine vert_interp_mixrat !------------------------------------------------------------------------------ subroutine vert_interp( ncol, levsiz, pin, pmid, datain, dataout ) !-------------------------------------------------------------------------- diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 2370cd416d..09df1a4214 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -2274,6 +2274,9 @@ 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 call ssatcontrail_d0(state, pbuf, ztodt, ptend, tend) call physics_update(state, ptend, ztodt, tend) diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index 2ca7096df6..3ba89bb04e 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -11,10 +11,11 @@ module ssatcontrail ! use phys_buffer, only: pbuf_size_max, pbuf_fld, pbuf_old_tim_idx, pbuf_get_fld_idx, pbuf_times use physics_buffer, only : pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_set_field, pbuf_old_tim_idx use constituents, only: cnst_get_ind, pcnst - use aircraft_emit, only: aircraft_cnt, spc_name_list use geopotential, only: geopotential_dse use phys_grid, only : get_wght_all_p use time_manager, only: get_curr_date + use aircraf_emit, only: get_aircraft + use wv_saturation, only: qsat_water, qsat_ice implicit none private @@ -39,8 +40,8 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) 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 - real(r8) :: w, esiT, ws, RH, ei + 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) @@ -52,6 +53,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) 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 @@ -61,16 +63,19 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) real(r8) :: dcld(pcols,pver) real(r8) :: ac_Q, ac_Q1, ac_Q2 real(r8) :: RHcts(pcols,pver) - integer :: itype logical :: lq(pcnst) integer :: yr, mon, day, ncsec real(r8) :: curr_factor - + integer :: aircraft_cnt + character(len=16) :: spc_name_list(30) + ! ICIWC = ICIWC0/rhodair has_aircraft_H2O = .false. + has_aircraft_distance = .false. + call get_aircraft(aircraft_cnt, spc_name_list) !----------------------------------------------------------------------------------------- ! check if ac_H2O in namelist ! if not, then bypass this subroutine @@ -80,10 +85,13 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) 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 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 @@ -153,10 +161,8 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) T_contr = T_contr + 273.15_r8 ! convert to Kelvin ! compute saturation pressure - itype = 0 - call gffgch(T_contr,eslTc,itype) - itype = 0 - call gffgch(state1%t(i,k),eslT,itype) + 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 @@ -164,8 +170,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) if(RH_contr.lt.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 - itype = 1 - call gffgch(state1%t(i,k),esiT,itype) ! saturation water vapor with respect to ice + 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) @@ -232,111 +237,6 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) enddo - call physics_update(state1, ptend_loc, dtime, tend) - end subroutine ssatcontrail_d0 -subroutine gffgch(t ,es ,itype ) - use shr_kind_mod, only: r8 => shr_kind_r8 - use physconst, only: tmelt - use abortutils, only: endrun - use cam_logfile, only: iulog - - implicit none -!------------------------------Arguments-------------------------------- -! -! Input arguments -! - real(r8), intent(in) :: t ! Temperature -! -! Output arguments -! - integer, intent(inout) :: itype ! Flag for ice phase and associated transition - - real(r8), intent(out) :: es ! Saturation vapor pressure -! -!---------------------------Local variables----------------------------- - real(r8) e1 ! Intermediate scratch variable for es over water - real(r8) e2 ! Intermediate scratch variable for es over water - real(r8) eswtr ! Saturation vapor pressure over water - real(r8) f ! Intermediate scratch variable for es over water - real(r8) f1 ! Intermediate scratch variable for es over water - real(r8) f2 ! Intermediate scratch variable for es over water - real(r8) f3 ! Intermediate scratch variable for es over water - real(r8) f4 ! Intermediate scratch variable for es over water - real(r8) f5 ! Intermediate scratch variable for es over water - real(r8) ps ! Reference pressure (mb) - real(r8) t0 ! Reference temperature (freezing point of water) - real(r8) term1 ! Intermediate scratch variable for es over ice - real(r8) term2 ! Intermediate scratch variable for es over ice - real(r8) term3 ! Intermediate scratch variable for es over ice - real(r8) tr ! Transition range for es over water to es over ice - real(r8) ts ! Reference temperature (boiling point of water) - real(r8) weight ! Intermediate scratch variable for es transition - integer itypo ! Intermediate scratch variable for holding itype -! -!----------------------------------------------------------------------- -! -! Check on whether there is to be a transition region for es -! - if (itype < 0) then - tr = abs(real(itype,r8)) - itypo = itype - itype = 1 - else - tr = 0.0_r8 - itypo = itype - end if - if (tr > 40.0_r8) then - write(iulog,900) tr - call endrun ('GFFGCH') ! Abnormal termination - end if -! - if(t < (tmelt - tr) .and. itype == 1) go to 10 -! -! Water -! - ps = 1013.246_r8 - ts = 373.16_r8 - e1 = 11.344_r8*(1.0_r8 - t/ts) - e2 = -3.49149_r8*(ts/t - 1.0_r8) - f1 = -7.90298_r8*(ts/t - 1.0_r8) - f2 = 5.02808_r8*log10(ts/t) - f3 = -1.3816_r8*(10.0_r8**e1 - 1.0_r8)/10000000.0_r8 - f4 = 8.1328_r8*(10.0_r8**e2 - 1.0_r8)/1000.0_r8 - f5 = log10(ps) - f = f1 + f2 + f3 + f4 + f5 - es = (10.0_r8**f)*100.0_r8 - eswtr = es -! - if(t >= tmelt .or. itype == 0) go to 20 -! -! Ice -! -10 continue - t0 = tmelt - term1 = 2.01889049_r8/(t0/t) - term2 = 3.56654_r8*log(t0/t) - term3 = 20.947031_r8*(t0/t) - es = 575.185606e10_r8*exp(-(term1 + term2 + term3)) -! - if (t < (tmelt - tr)) go to 20 -! -! Weighted transition between water and ice -! - weight = min((tmelt - t)/tr,1.0_r8) - es = weight*es + (1.0_r8 - weight)*eswtr -! -20 continue - itype = itypo - return -! -900 format('GFFGCH: FATAL ERROR ******************************',/, & - 'TRANSITION RANGE FOR WATER TO ICE SATURATION VAPOR', & - ' PRESSURE, TR, EXCEEDS MAXIMUM ALLOWABLE VALUE OF', & - ' 40.0 DEGREES C',/, ' TR = ',f7.2) -! -end subroutine gffgch - - end module ssatcontrail From 7f88f7126ccf4928a5a8cbd42ad30e17d80012db Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Thu, 14 Jan 2021 00:14:55 -0700 Subject: [PATCH 04/24] revised contrail code 1/14/2021 --- .../utils/horizontal_interpolate.F90 | 23 ++++++++++--------- src/chemistry/utils/tracer_data.F90 | 12 ++++++---- src/physics/cam/ssatcontrail.F90 | 2 -- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/chemistry/utils/horizontal_interpolate.F90 b/src/chemistry/utils/horizontal_interpolate.F90 index c2d3212b64..2d2a72930c 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,ii) + subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) !------------------------------------------------------------------------------------------------------------ ! 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 @@ -27,7 +27,8 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! !------------------------------------------------------------------------------------------------------------ implicit none - integer, intent(in) :: im1, jm1, im2, jm2, ii + integer, intent(in) :: im1, jm1, im2, jm2 + logical, intent(in) :: itype !.true. = flight distance only, .false. = all mixing ratios real(r8), intent(in) :: lon0(im1), lat0(jm1) real(r8), intent(out) :: weight_x(im2,im1), weight_y(jm2,jm1) @@ -116,7 +117,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |-------------------| ! |---------------------------------| ! x2_west x2_east - if(ii.eq.0) then + if(itype.eq.0) then weight_x(i2,i1) = 1.0_r8 else weight_x(i2,i1) = (x1_east-x1_west)/(x2_east-x2_west) @@ -127,7 +128,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |--------------------------------| ! |---------------------------------| ! x2_west x2_east - if(ii.eq.0) then + if(itype.eq.0) 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) @@ -138,7 +139,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |--------------------------------| ! |---------------------------------| ! x2_west x2_east - if(ii.eq.0) then + if(itype.eq.0) 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) @@ -163,7 +164,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |-------------------------| ! |----------------|......................| ! slon2(im2) slon2(im2+1) slon2(2) (note: slon2(im2+1) = slon2(1)) - if(ii.eq.0) then + if(itype.eq.0) 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)) @@ -176,7 +177,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |-------------------------|.............................| ! |-------------------------------| ! slon2(im2) slon2(im2+1) <--- end point - if(ii.eq.0) then + if(itype.eq.0) 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)) @@ -207,7 +208,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |-------------------| ! |---------------------------------| ! y2_south y2_north - if(ii.eq.0) then + if(itype.eq.0) then weight_y(j2,j1) = 1.0_r8 else weight_y(j2,j1) = gw1(j1)/gw2(j2) @@ -218,7 +219,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |--------------------------------| ! |---------------------------------| ! y2_south y2_north - if(ii.eq.0) then + if(itype.eq.0) 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) @@ -229,7 +230,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |--------------------------------| ! |---------------------------------| ! y2_south y2_north - if(ii.eq.0) then + if(itype.eq.0) 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) @@ -240,7 +241,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,ii) ! |--------------------------------| ! |---------------------| ! y2_south y2_north - if(ii.eq.0) then + if(itype.eq.0) then weight_y(j2,j1) = 1._r8 else weight_y(j2,j1) = gw1(j1)/gw2(j2) diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index e72429c253..1c0c6650ba 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -628,7 +628,8 @@ 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)) + if(aircraft_cnt.gt.0) then + allocate(file%weight_x(plon,file%nlon)) allocate(file%weight_y(plat,file%nlat)) allocate(file%count_x(plon)) allocate(file%count_y(plat)) @@ -656,7 +657,8 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & 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,1) + call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats, & + plon,plat,file%weight_x,file%weight_y,.false.) do i2=1,plon file%count_x(i2) = 0 @@ -678,7 +680,8 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & enddo enddo - call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats,plon,plat,file%weight0_x,file%weight0_y,0) + call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats,& + plon,plat,file%weight0_x,file%weight0_y,.true.) do i2=1,plon file%count0_x(i2) = 0 @@ -700,7 +703,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & 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) @@ -715,6 +718,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & call mpibcast(file%index0_x, plon*file%nlon, mpiint , 0, mpicom) call mpibcast(file%index0_y, plat*file%nlat, mpiint , 0, mpicom) #endif + endif endif end subroutine trcdata_init diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index 3ba89bb04e..2cdccfaf45 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -8,7 +8,6 @@ module ssatcontrail use physics_types, only: physics_ptend_sum, physics_update use physics_types, only: physics_state_copy, physics_ptend_init use physconst, only: cpair,mwdry,mwh2o, gravit, zvir, rair, pi, rearth -! use phys_buffer, only: pbuf_size_max, pbuf_fld, pbuf_old_tim_idx, pbuf_get_fld_idx, pbuf_times use physics_buffer, only : pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_set_field, pbuf_old_tim_idx use constituents, only: cnst_get_ind, pcnst use geopotential, only: geopotential_dse @@ -24,7 +23,6 @@ module ssatcontrail public ssatcontrail_d0 ! Private data -! real(r8), parameter :: ICIWC0 = 2.0e-6 ! ICIWC = 2 mg/m3, converted to kg/kg here by dividing rho_air real(r8), parameter :: rhoi = 500.0_r8 ! density of ice (500 kg/m3) real(r8), parameter :: radius = 3.75e-6 ! diameter of ice particle = 7.5 microns From 78065ab339e14c055ac19bb8215a3c3bafc5f038 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Thu, 14 Jan 2021 21:53:24 -0700 Subject: [PATCH 05/24] revision 1/15/2021 --- .../utils/horizontal_interpolate.F90 | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/chemistry/utils/horizontal_interpolate.F90 b/src/chemistry/utils/horizontal_interpolate.F90 index 2d2a72930c..1c854440de 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,itype) + 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,7 +28,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) !------------------------------------------------------------------------------------------------------------ implicit none integer, intent(in) :: im1, jm1, im2, jm2 - logical, intent(in) :: itype !.true. = flight distance only, .false. = all mixing ratios + 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) @@ -117,7 +117,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |-------------------| ! |---------------------------------| ! x2_west x2_east - if(itype.eq.0) then + 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) @@ -128,7 +128,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |--------------------------------| ! |---------------------------------| ! x2_west x2_east - if(itype.eq.0) then + 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) @@ -139,7 +139,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |--------------------------------| ! |---------------------------------| ! x2_west x2_east - if(itype.eq.0) then + 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) @@ -164,7 +164,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |-------------------------| ! |----------------|......................| ! slon2(im2) slon2(im2+1) slon2(2) (note: slon2(im2+1) = slon2(1)) - if(itype.eq.0) then + 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)) @@ -177,7 +177,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |-------------------------|.............................| ! |-------------------------------| ! slon2(im2) slon2(im2+1) <--- end point - if(itype.eq.0) then + 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)) @@ -208,7 +208,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |-------------------| ! |---------------------------------| ! y2_south y2_north - if(itype.eq.0) then + if(use_flight_distance) then weight_y(j2,j1) = 1.0_r8 else weight_y(j2,j1) = gw1(j1)/gw2(j2) @@ -219,7 +219,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |--------------------------------| ! |---------------------------------| ! y2_south y2_north - if(itype.eq.0) 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) @@ -230,7 +230,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |--------------------------------| ! |---------------------------------| ! y2_south y2_north - if(itype.eq.0) 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) @@ -241,7 +241,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,itype) ! |--------------------------------| ! |---------------------| ! y2_south y2_north - if(itype.eq.0) then + if(use_flight_distance) then weight_y(j2,j1) = 1._r8 else weight_y(j2,j1) = gw1(j1)/gw2(j2) From b27dc544437fb34e88c4e0ca3d6dc6231399af42 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Thu, 14 Jan 2021 22:02:45 -0700 Subject: [PATCH 06/24] revision 1/15/2021 v2 --- src/chemistry/utils/tracer_data.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 1c0c6650ba..4ec10a3848 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -1902,11 +1902,11 @@ subroutine interpolate_trcdata( state, flds, file, pbuf2d ) if( file%dist ) then call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, & datain, data_out(:,:), & - file%p0,ps,file%hyai,file%hybi,0) + file%p0,ps,file%hyai,file%hybi,.true.) else call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, & datain, data_out(:,:), & - file%p0,ps,file%hyai,file%hybi,1) + file%p0,ps,file%hyai,file%hybi,.false.) endif else call vert_interp(ncol, file%nlev, pin, state(c)%pmid, datain, data_out(:,:) ) @@ -2378,7 +2378,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, ii) + subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, hybi, use_flight_distance) implicit none @@ -2388,7 +2388,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 - integer, intent(in) :: ii ! 0: rebin only, otherwise mixing ratio + logical, intent(in) :: use_flight_distance ! .true. = flight distance, .flase. = mixing ratio real(r8), intent(out) :: trg(pcols,ntrg) ! target array real(r8) :: ps(pcols), p0, hyai(nsrc+1), hybi(nsrc+1) @@ -2424,14 +2424,14 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, top = trg_x(n,i) do j = sil,1,-1 if( top.lt.src_x(j) ) then - if(ii.eq.0) then + if(use_flight_distance) then y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) else y = y+(bot-src_x(j))*src(n,j) endif bot = src_x(j) else - if(ii.eq.0) then + if(use_flight_distance) then y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) else y = y+(bot-top)*src(n,j) @@ -2451,14 +2451,14 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, y = 0.0_r8 do j=nsrc,1,-1 if( top.lt.src_x(j) ) then - if(ii.eq.0) then + if(use_flight_distance) then y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) else y = y+(bot-src_x(j))*src(n,j) endif bot = src_x(j) else - if(ii.eq.0) then + if(use_flight_distance) then y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) else y = y+(bot-top)*src(n,j) @@ -2470,7 +2470,7 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, endif ! turn mass into mixing ratio - if(ii.ne.0) then + if(.not. use_flight_distance) then do i=1,ntrg trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i)) enddo From 38434609662a5e3483b2aca01240cda1ae0c0b44 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Mon, 18 Jan 2021 01:09:14 -0700 Subject: [PATCH 07/24] revision 1/18/2021 --- src/chemistry/utils/tracer_data.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 4ec10a3848..ee873351ee 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -1685,8 +1685,9 @@ subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order) 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%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) + 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 @@ -1694,8 +1695,9 @@ subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order) 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) + 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') From 381766f3cf945e9aea752f67f8a12a34e72c1917 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Sun, 24 Jan 2021 21:13:10 -0700 Subject: [PATCH 08/24] revised 1/25/2021 --- src/chemistry/utils/aircraft_emit.F90 | 2 +- src/chemistry/utils/tracer_data.F90 | 10 ++++++---- src/physics/cam/ssatcontrail.F90 | 24 ++++++++++++------------ 3 files changed, 19 insertions(+), 17 deletions(-) diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index ea05fa58d8..05e3298e6f 100644 --- a/src/chemistry/utils/aircraft_emit.F90 +++ b/src/chemistry/utils/aircraft_emit.F90 @@ -78,7 +78,7 @@ subroutine get_aircraft(cnt, spc_name_list_out) integer :: i cnt = aircraft_cnt - if( cnt.gt. 0 ) then + if( cnt>0 ) then do i=1,cnt spc_name_list_out(i) = spc_name_list(i) end do diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index ee873351ee..b9b288f98d 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -597,7 +597,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & call get_horiz_grid_d(plon, clon_d_out=lam) call get_aircraft(aircraft_cnt) - if(aircraft_cnt.gt.0) then + if(aircraft_cnt>0) then 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)) endif @@ -628,7 +628,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & deallocate(phi,lam) ! weight_x & weight_y are weighting function for x & y interpolation - if(aircraft_cnt.gt.0) then + if(aircraft_cnt>0) then allocate(file%weight_x(plon,file%nlon)) allocate(file%weight_y(plat,file%nlat)) allocate(file%count_x(plon)) @@ -703,8 +703,10 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & enddo enddo endif - + endif + #if ( defined SPMD) + if(aircraft_cnt>0) then 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) @@ -717,8 +719,8 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & call mpibcast(file%count0_y, plat, mpiint , 0, mpicom) call mpibcast(file%index0_x, plon*file%nlon, mpiint , 0, mpicom) call mpibcast(file%index0_y, plat*file%nlat, mpiint , 0, mpicom) + endif #endif - endif endif end subroutine trcdata_init diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index 2cdccfaf45..d6c0793bb3 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -1,5 +1,7 @@ module ssatcontrail - +! contrail parameterization +! see Chen et al., 2012: Global contrail coverage simulated +! by CAM5 with the inventory of 2006 global aircraft emissions, JAMES use shr_kind_mod, only: r8 => shr_kind_r8 use ppgrid, only: pcols, pver use cam_history, only: outfld @@ -68,8 +70,6 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) integer :: aircraft_cnt character(len=16) :: spc_name_list(30) -! ICIWC = ICIWC0/rhodair - has_aircraft_H2O = .false. has_aircraft_distance = .false. @@ -91,6 +91,15 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) if(.not. has_aircraft_H2O) return if(.not. has_aircraft_distance) return +!------------------------------------------------------------------------------------------ + lq(:) = .FALSE. + lq(1) = .TRUE. + lq(ixcldice) = .TRUE. + lq(ixnumice) = .TRUE. + + call physics_ptend_init(ptend_loc, state1%psetcols,'ssatcontrail',ls=.true.,lq = lq) +!----------------------------------------------------------------------------------------- + particle_mass = 4._r8/3._r8*pi*rhoi*radius**3 ! mass of ice particle rog = rair/gravit ! Rd/Cp @@ -124,15 +133,6 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) call cnst_get_ind('NUMICE', ixnumice, abort=.false.) call cnst_get_ind('NUMLIQ', ixnumliq, abort=.false.) -!------------------------------------------------------------------------------------------ - lq(:) = .FALSE. - lq(1) = .TRUE. - lq(ixcldice) = .TRUE. - lq(ixnumice) = .TRUE. - - call physics_ptend_init(ptend_loc, state1%psetcols,'ssatcontrail',ls=.true.,lq = lq) -!----------------------------------------------------------------------------------------- - ! adjust h2o to volume mixing ratio (mass adjustment and conversion from g/kg to kg/kg) Ma = mwdry Mh2o = mwh2o From 43c2be5aaca3881c74bd643e37b4defc49e77903 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Thu, 11 Feb 2021 00:18:07 -0700 Subject: [PATCH 09/24] revised 2/11/2021 --- src/chemistry/utils/tracer_data.F90 | 38 ++++++++++++++--------------- src/physics/cam/ssatcontrail.F90 | 8 +++--- 2 files changed, 22 insertions(+), 24 deletions(-) diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index b9b288f98d..10fa46a634 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -171,7 +171,6 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & #if ( defined SPMD ) use mpishorthand, only: mpicom, mpir8, mpiint #endif - use aircraf_emit, only: get_aircraft implicit none @@ -204,7 +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 :: aircrat_cnt + integer :: aircraft_cnt call specify_fields( specifier, flds ) @@ -596,11 +595,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) - call get_aircraft(aircraft_cnt) - if(aircraft_cnt>0) then 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)) - endif lon_global_grid_ndx=-huge(1) lat_global_grid_ndx=-huge(1) @@ -628,7 +624,6 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & deallocate(phi,lam) ! weight_x & weight_y are weighting function for x & y interpolation - if(aircraft_cnt>0) then allocate(file%weight_x(plon,file%nlon)) allocate(file%weight_y(plat,file%nlat)) allocate(file%count_x(plon)) @@ -642,18 +637,20 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & file%index_x(:,:) = 0 file%index_y(:,:) = 0 - allocate(file%weight0_x(plon,file%nlon)) - allocate(file%weight0_y(plat,file%nlat)) - allocate(file%count0_x(plon)) - allocate(file%count0_y(plat)) - allocate(file%index0_x(plon,file%nlon)) - allocate(file%index0_y(plat,file%nlat)) - 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 + if( file%dist ) then + allocate(file%weight0_x(plon,file%nlon)) + allocate(file%weight0_y(plat,file%nlat)) + allocate(file%count0_x(plon)) + allocate(file%count0_y(plat)) + allocate(file%index0_x(plon,file%nlon)) + allocate(file%index0_y(plat,file%nlat)) + 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 @@ -680,6 +677,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & enddo enddo + if( file%dist ) then call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats,& plon,plat,file%weight0_x,file%weight0_y,.true.) @@ -702,17 +700,17 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & endif enddo enddo + endif endif - endif #if ( defined SPMD) - if(aircraft_cnt>0) then 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) + if( file%dist ) then call mpibcast(file%weight0_x, plon*file%nlon, mpir8 , 0, mpicom) call mpibcast(file%weight0_y, plat*file%nlat, mpir8 , 0, mpicom) call mpibcast(file%count0_x, plon, mpiint , 0, mpicom) diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index d6c0793bb3..5cc242761b 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -2,6 +2,7 @@ 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 cam_history, only: outfld @@ -15,8 +16,8 @@ module ssatcontrail use geopotential, only: geopotential_dse use phys_grid, only : get_wght_all_p use time_manager, only: get_curr_date - use aircraf_emit, only: get_aircraft use wv_saturation, only: qsat_water, qsat_ice + use aircraft_emit, only: get_aircraft implicit none private @@ -88,9 +89,6 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) endif enddo endif - if(.not. has_aircraft_H2O) return - if(.not. has_aircraft_distance) return - !------------------------------------------------------------------------------------------ lq(:) = .FALSE. lq(1) = .TRUE. @@ -99,6 +97,8 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) 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 From f61f53a659103f72b341d6f760c43d7dd94639d3 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 1 Mar 2021 17:17:48 -0700 Subject: [PATCH 10/24] Add error handling and update MPI --- src/chemistry/utils/tracer_data.F90 | 120 +++++++++++++++++++++------- 1 file changed, 89 insertions(+), 31 deletions(-) diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 8beac39415..47236a1d70 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -168,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 @@ -186,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 @@ -624,12 +624,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 @@ -638,12 +662,36 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & file%index_y(:,:) = 0 if( file%dist ) then - allocate(file%weight0_x(plon,file%nlon)) - allocate(file%weight0_y(plat,file%nlat)) - allocate(file%count0_x(plon)) - allocate(file%count0_y(plat)) - allocate(file%index0_x(plon,file%nlon)) - allocate(file%index0_y(plat,file%nlat)) + 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 @@ -703,22 +751,32 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & endif 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) - if( file%dist ) then - call mpibcast(file%weight0_x, plon*file%nlon, mpir8 , 0, mpicom) - call mpibcast(file%weight0_y, plat*file%nlat, mpir8 , 0, mpicom) - call mpibcast(file%count0_x, plon, mpiint , 0, mpicom) - call mpibcast(file%count0_y, plat, mpiint , 0, mpicom) - call mpibcast(file%index0_x, plon*file%nlon, mpiint , 0, mpicom) - call mpibcast(file%index0_y, plat*file%nlat, mpiint , 0, mpicom) - endif -#endif + call mpi_bcast(file%weight_x, 1, mpi_real8 , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight_x") + call mpi_bcast(file%weight_y, 1, mpi_real8 , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight_y") + call mpi_bcast(file%count_x, 1, mpi_integer , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count_x") + call mpi_bcast(file%count_y, 1, mpi_integer , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count_y") + call mpi_bcast(file%index_x, 1, mpi_integer , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index_x") + call mpi_bcast(file%index_y, 1, 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, 1, mpi_real8 , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight0_x") + call mpi_bcast(file%weight0_y, 1, mpi_real8 , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%weight0_y") + call mpi_bcast(file%count0_x, 1, mpi_integer , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count0_x") + call mpi_bcast(file%count0_y, 1, mpi_integer , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%count0_y") + call mpi_bcast(file%index0_x, 1, mpi_integer , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index0_x") + call mpi_bcast(file%index0_y, 1, mpi_integer , mstrid, mpicom,ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: file%index0_y") + endif endif end subroutine trcdata_init From fa7b53de66a5d3967daf960df29ed637a94cdda9 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Wed, 3 Mar 2021 15:19:33 -0700 Subject: [PATCH 11/24] Fix errors from regression testing --- src/chemistry/utils/aircraft_emit.F90 | 6 +++--- src/physics/cam/ssatcontrail.F90 | 22 +++++++++++----------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index 05e3298e6f..f295066d53 100644 --- a/src/chemistry/utils/aircraft_emit.F90 +++ b/src/chemistry/utils/aircraft_emit.F90 @@ -42,9 +42,9 @@ module aircraft_emit type(forcing_air),allocatable :: forcings_air(:) 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 '/) + 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 diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index 5cc242761b..15d31f7582 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -27,7 +27,7 @@ module ssatcontrail ! Private data real(r8), parameter :: rhoi = 500.0_r8 ! density of ice (500 kg/m3) - real(r8), parameter :: radius = 3.75e-6 ! diameter of ice particle = 7.5 microns + real(r8), parameter :: radius = 3.75e-6_r8 ! diameter of ice particle = 7.5 microns contains @@ -74,6 +74,14 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) 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 @@ -125,14 +133,6 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) call pbuf_get_field(pbuf,ifld,ac_SLANT_DIST) - ! 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.) - ! adjust h2o to volume mixing ratio (mass adjustment and conversion from g/kg to kg/kg) Ma = mwdry Mh2o = mwh2o @@ -142,7 +142,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) epsi = Mh2o/Ma ei = 1.21_r8 ! water vapor emiision index (g) h2o per kg fuel (Schumann 96)? - Q = 43.e6 ! specific combustion heat Schummann 1996, Q = 43 MJ/kg + Q = 43.e6_r8 ! specific combustion heat Schummann 1996, Q = 43 MJ/kg eta = 0.3_r8 ! propulsion effieciency (Ponater 2002) ratio(i,k) = 0._r8 @@ -180,7 +180,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) ICIWC0(i,k) = exp(6.97_r8+0.103_r8*(state1%t(i,k)-273.16_r8)) ! in mg/m3 rho = p/(287._r8*state1%t(i,k)) - ICIWC = ICIWC0(i,k)/rho*1.0e-6 + ICIWC = ICIWC0(i,k)/rho*1.0e-6_r8 ! persistent contrail condition From 70e0d73924a412d9d54308e5b75d7e41ba5c1d7d Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Thu, 4 Mar 2021 23:55:18 -0700 Subject: [PATCH 12/24] revised 3/5/2021 --- src/chemistry/utils/aircraft_emit.F90 | 8 +-- src/chemistry/utils/tracer_data.F90 | 76 +++++++++++++-------------- src/physics/cam/physpkg.F90 | 1 + src/physics/cam/ssatcontrail.F90 | 10 ++-- 4 files changed, 48 insertions(+), 47 deletions(-) diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index f295066d53..876edb9147 100644 --- a/src/chemistry/utils/aircraft_emit.F90 +++ b/src/chemistry/utils/aircraft_emit.F90 @@ -68,7 +68,7 @@ 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) - integer :: dist(N_AERO) + logical :: dist(N_AERO) contains @@ -106,7 +106,7 @@ subroutine aircraft_emit_register() !------------------------------------------------------------------ ! Return if air_specifier is blank (no aircraft data to process) !------------------------------------------------------------------ - dist(:) = 0 + dist(:) = .false. aircraft_cnt = 0 if (air_specifier(1) == "") return @@ -126,7 +126,7 @@ 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) = 1 + 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) @@ -209,7 +209,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 - if( dist(m) == 1 ) forcings_air(m)%file%dist = .true. + 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 diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 47236a1d70..eaf0a9fab4 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -662,42 +662,42 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & 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 + 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 @@ -2451,7 +2451,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, .flase. = mixing ratio + 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) @@ -2537,7 +2537,7 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, do i=1,ntrg trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i)) enddo - endif + endif enddo diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 95c66b09fe..b6bd7780bb 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -2277,6 +2277,7 @@ subroutine tphysbc (ztodt, state, & ! 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, tend) call physics_update(state, ptend, ztodt, tend) diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index 15d31f7582..3b46d254e0 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -10,7 +10,7 @@ module ssatcontrail use physics_types, only: physics_state, physics_ptend, physics_tend use physics_types, only: physics_ptend_sum, physics_update use physics_types, only: physics_state_copy, physics_ptend_init - use physconst, only: cpair,mwdry,mwh2o, gravit, zvir, rair, pi, rearth + 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_set_field, pbuf_old_tim_idx use constituents, only: cnst_get_ind, pcnst use geopotential, only: geopotential_dse @@ -35,7 +35,7 @@ module ssatcontrail subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) implicit none - type(physics_state), intent(inout) :: state1 + type(physics_state), intent(in) :: state1 type(physics_ptend), intent(inout) :: ptend_loc type(physics_tend) :: tend type(physics_buffer_desc), pointer :: pbuf(:) @@ -156,7 +156,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) 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 + 273.15_r8 ! convert to Kelvin + T_contr = T_contr + tmelt ! convert to Kelvin ! compute saturation pressure call qsat_water(T_contr, p, eslTc, qslTc) @@ -178,8 +178,8 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) ! Schumann, 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)-273.16_r8)) ! in mg/m3 - rho = p/(287._r8*state1%t(i,k)) + 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 From 04335eae253925f7aaef6d4e05966c6ab11eb213 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 8 Mar 2021 08:54:18 -0700 Subject: [PATCH 13/24] Additional changes based on reviewer comments --- .../utils/horizontal_interpolate.F90 | 36 ++--- src/chemistry/utils/tracer_data.F90 | 126 +++++++++--------- 2 files changed, 81 insertions(+), 81 deletions(-) diff --git a/src/chemistry/utils/horizontal_interpolate.F90 b/src/chemistry/utils/horizontal_interpolate.F90 index 1c854440de..00f3a6206e 100644 --- a/src/chemistry/utils/horizontal_interpolate.F90 +++ b/src/chemistry/utils/horizontal_interpolate.F90 @@ -164,11 +164,11 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight ! |-------------------------| ! |----------------|......................| ! slon2(im2) slon2(im2+1) slon2(2) (note: slon2(im2+1) = 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 + 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 @@ -177,11 +177,11 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight ! |-------------------------|.............................| ! |-------------------------------| ! slon2(im2) slon2(im2+1) <--- end point - 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 + 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 @@ -209,9 +209,9 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight ! |---------------------------------| ! y2_south y2_north if(use_flight_distance) then - weight_y(j2,j1) = 1.0_r8 + weight_y(j2,j1) = 1.0_r8 else - weight_y(j2,j1) = gw1(j1)/gw2(j2) + weight_y(j2,j1) = gw1(j1)/gw2(j2) endif elseif ( (y1_south.ge.y2_south).and.(y1_south.lt.y2_north) ) then ! case 2: @@ -220,9 +220,9 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight ! |---------------------------------| ! y2_south y2_north if(use_flight_distance) then - weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south) + 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) + weight_y(j2,j1) = (y2_north-y1_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2) endif elseif ( (y1_north.gt.y2_south).and.(y1_north.le.y2_north) ) then ! case 3: @@ -231,9 +231,9 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight ! |---------------------------------| ! y2_south y2_north if(use_flight_distance) then - weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south) + 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) + weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2) endif elseif ( (y1_north.gt.y2_north).and.(y1_south.lt.y2_south) ) then ! case 4: @@ -242,9 +242,9 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight ! |---------------------| ! y2_south y2_north if(use_flight_distance) then - weight_y(j2,j1) = 1._r8 + weight_y(j2,j1) = 1._r8 else - weight_y(j2,j1) = gw1(j1)/gw2(j2) + weight_y(j2,j1) = gw1(j1)/gw2(j2) endif endif diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index eaf0a9fab4..37eb351d77 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -2467,77 +2467,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,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, 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 if( tl.gt.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.lt.src_x(j) ) then + if(use_flight_distance) then + y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) + else + y = y+(bot-src_x(j))*src(n,j) + endif + bot = src_x(j) + else + if(use_flight_distance) then + y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) + else + y = y+(bot-top)*src(n,j) + endif + exit + endif + enddo + trg(n,i) = y + else + trg(n,i) = 0.0_r8 + end if + end do + + if( trg_x(n,ntrg+1).lt.src_x(nsrc+1) ) then + top = trg_x(n,ntrg+1) + bot = src_x(nsrc+1) y = 0.0_r8 - bot = min(tl,src_x(nsrc+1)) - top = trg_x(n,i) - do j = sil,1,-1 - if( top.lt.src_x(j) ) then - if(use_flight_distance) then - y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) - else - y = y+(bot-src_x(j))*src(n,j) - endif - bot = src_x(j) - else - if(use_flight_distance) then - y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) - else - y = y+(bot-top)*src(n,j) - endif - exit - endif + do j=nsrc,1,-1 + if( top.lt.src_x(j) ) then + if(use_flight_distance) then + y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) + else + y = y+(bot-src_x(j))*src(n,j) + endif + bot = src_x(j) + else + if(use_flight_distance) then + y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) + else + y = y+(bot-top)*src(n,j) + endif + exit + endif enddo - trg(n,i) = y - else - trg(n,i) = 0.0_r8 - end if - end do - - if( trg_x(n,ntrg+1).lt.src_x(nsrc+1) ) then - top = trg_x(n,ntrg+1) - bot = src_x(nsrc+1) - y = 0.0_r8 - do j=nsrc,1,-1 - if( top.lt.src_x(j) ) then - if(use_flight_distance) then - y = y+(bot-src_x(j))*src(n,j)/(src_x(j+1)-src_x(j)) - else - y = y+(bot-src_x(j))*src(n,j) - endif - bot = src_x(j) - else - if(use_flight_distance) then - y = y+(bot-top)*src(n,j)/(src_x(j+1)-src_x(j)) - else - y = y+(bot-top)*src(n,j) + trg(n,ntrg) = trg(n,ntrg)+y endif - exit - endif - enddo - trg(n,ntrg) = trg(n,ntrg)+y - endif -! turn mass into mixing ratio - if(.not. use_flight_distance) then - do i=1,ntrg - trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i)) - enddo - endif + ! turn mass into mixing ratio + if(.not. use_flight_distance) then + do i=1,ntrg + trg(n,i) = trg(n,i)/(trg_x(n,i+1)-trg_x(n,i)) + enddo + endif enddo From 7202872d5a23656919c83cf78fcfdfa25e473ab5 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 8 Mar 2021 09:04:10 -0700 Subject: [PATCH 14/24] Add DOI --- src/physics/cam/ssatcontrail.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index 3b46d254e0..07f2291c01 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -138,7 +138,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) Mh2o = mwh2o ! contrail paramter (G = CF*p/epi) -! and Schumann 1996, reprinted by Ponater 2002, JGR (eq 6-8) +! and Schumann 1996, reprinted by Ponater 2002, JGR (eq 6-8) DOI: 10.1029/2011MS000105 epsi = Mh2o/Ma ei = 1.21_r8 ! water vapor emiision index (g) h2o per kg fuel (Schumann 96)? From 4d02c1c785fa1f03d0c877ed1fe01ed11364da02 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 8 Mar 2021 09:16:12 -0700 Subject: [PATCH 15/24] Add book citation --- src/physics/cam/ssatcontrail.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index 07f2291c01..91c19b143f 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -175,7 +175,8 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) RH = w/ws ! relative humidity with respect to ice if( RH.ge.1.0_r8 ) RHcts(i,k) = 1.0_r8 -! Schumann, 2002: IWC(g/m3) = exp(6.97+0.103*T(C))*1e-3 +! 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 From ae9186e1778644f8f202386cb3f4c60fee778d62 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Mon, 8 Mar 2021 19:20:18 -0700 Subject: [PATCH 16/24] revised 3/9/2021 --- src/chemistry/utils/aircraft_emit.F90 | 5 ++--- src/physics/cam/ssatcontrail.F90 | 16 ++++++++-------- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index 876edb9147..56eb15e412 100644 --- a/src/chemistry/utils/aircraft_emit.F90 +++ b/src/chemistry/utils/aircraft_emit.F90 @@ -48,9 +48,8 @@ module aircraft_emit real(r8), parameter :: molmass(N_AERO) = 1._r8 - logical :: advective_tracer(N_AERO) = (/.false., .false., .false., .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','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 diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index 91c19b143f..ebb579de01 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -11,13 +11,13 @@ module ssatcontrail use physics_types, only: physics_ptend_sum, physics_update use physics_types, only: physics_state_copy, 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_set_field, pbuf_old_tim_idx - use constituents, only: cnst_get_ind, pcnst - use geopotential, only: geopotential_dse - use phys_grid, only : get_wght_all_p - use time_manager, only: get_curr_date - use wv_saturation, only: qsat_water, qsat_ice - use aircraft_emit, only: get_aircraft + use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_set_field, pbuf_old_tim_idx + use constituents, only: cnst_get_ind, pcnst + use geopotential, only: geopotential_dse + use phys_grid, only: get_wght_all_p + use time_manager, only: get_curr_date + use wv_saturation, only: qsat_water, qsat_ice + use aircraft_emit, only: get_aircraft implicit none private @@ -138,7 +138,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) Mh2o = mwh2o ! contrail paramter (G = CF*p/epi) -! and Schumann 1996, reprinted by Ponater 2002, JGR (eq 6-8) DOI: 10.1029/2011MS000105 +! 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 emiision index (g) h2o per kg fuel (Schumann 96)? From 822e2638fcf08756426d381b85ec5c27eaa4a9a1 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Mon, 8 Mar 2021 19:29:03 -0700 Subject: [PATCH 17/24] revised 3/9/2021 v2 --- src/chemistry/utils/horizontal_interpolate.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/chemistry/utils/horizontal_interpolate.F90 b/src/chemistry/utils/horizontal_interpolate.F90 index 00f3a6206e..de0602213a 100644 --- a/src/chemistry/utils/horizontal_interpolate.F90 +++ b/src/chemistry/utils/horizontal_interpolate.F90 @@ -235,7 +235,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight else weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2) endif - elseif ( (y1_north.gt.y2_north).and.(y1_south.lt.y2_south) ) then + elseif ( (y1_north.gt.y2_north).and.(y1_south.lt.y2_south) ) then ! case 4: ! y1_south y1_north ! |--------------------------------| From 5ee82946462db5e42b55607957512ef4e32f8c08 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Mon, 8 Mar 2021 22:07:19 -0700 Subject: [PATCH 18/24] revised 3/9/2021 v3 --- src/chemistry/utils/aircraft_emit.F90 | 2 ++ src/chemistry/utils/tracer_data.F90 | 25 ++++++++++--------------- src/physics/cam/physpkg.F90 | 2 +- src/physics/cam/ssatcontrail.F90 | 19 +++++++------------ 4 files changed, 20 insertions(+), 28 deletions(-) diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index 56eb15e412..f2a242aa63 100644 --- a/src/chemistry/utils/aircraft_emit.F90 +++ b/src/chemistry/utils/aircraft_emit.F90 @@ -76,6 +76,8 @@ subroutine get_aircraft(cnt, spc_name_list_out) 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 diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 37eb351d77..b1528dc811 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -703,7 +703,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & 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,.false.) + plon,plat,file%weight_x,file%weight_y,file%dist) do i2=1,plon file%count_x(i2) = 0 @@ -727,12 +727,12 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & if( file%dist ) then call xy_interp_init(file%nlon,file%nlat,file%lons,file%lats,& - plon,plat,file%weight0_x,file%weight0_y,.true.) + 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).gt.0.0_r8 ) then + 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 @@ -742,7 +742,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & do j2=1,plat file%count0_y(j2) = 0 do j1=1,file%nlat - if(file%weight0_y(j2,j1).gt.0.0_r8 ) then + 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 @@ -1962,14 +1962,9 @@ subroutine interpolate_trcdata( state, flds, file, pbuf2d ) if ( file%top_bndry ) then call vert_interp_ub(ncol, file%nlev, file%levs, datain(:ncol,:), data_out(:ncol,:) ) else if(file%conserve_column) then - if( file%dist ) then call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, & datain, data_out(:,:), & - file%p0,ps,file%hyai,file%hybi,.true.) - else - call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, & - datain, data_out(:,:), & - file%p0,ps,file%hyai,file%hybi,.false.) + file%p0,ps,file%hyai,file%hybi,file%dist) endif else call vert_interp(ncol, file%nlev, pin, state(c)%pmid, datain, data_out(:,:) ) @@ -2473,9 +2468,9 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, 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 + if( (tl>src_x(1)).and.(trg_x(n,i) shr_kind_r8 use ppgrid, only: pcols, pver - use cam_history, only: outfld - use cam_logfile, only: iulog - use physics_types, only: physics_state, physics_ptend, physics_tend - use physics_types, only: physics_ptend_sum, physics_update - use physics_types, only: physics_state_copy, physics_ptend_init + 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_set_field, pbuf_old_tim_idx use constituents, only: cnst_get_ind, pcnst @@ -32,12 +28,11 @@ module ssatcontrail contains - subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) + 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_tend) :: tend type(physics_buffer_desc), pointer :: pbuf(:) real(r8), intent(in) :: dtime ! time step !------------------------Local storage------------------------------------------------------ @@ -141,7 +136,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) ! 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 emiision index (g) h2o per kg fuel (Schumann 96)? + 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) @@ -164,8 +159,8 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) RH_contr = (G*(state1%t(i,k)-T_contr)+eslTc)/eslT ! RH_contr ranges between 0 and 1 - if(RH_contr.gt.1.0_r8) RH_contr = 1.0_r8 - if(RH_contr.lt.0.0_r8) RH_contr = 0.0_r8 + 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) @@ -173,7 +168,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) qs = ws/(1.0_r8+ws) RH = w/ws ! relative humidity with respect to ice - if( RH.ge.1.0_r8 ) RHcts(i,k) = 1.0_r8 + 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 @@ -185,7 +180,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc, tend) ! persistent contrail condition - if( (state1%t(i,k).lt.T_contr).and.(RH.gt.RH_contr).and.(RH.gt.1.0_r8).and.(ac_H2O(i,k).gt.0.0_r8) ) then + 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) From ee7201612674b81131aa3e86757c9b2438c3a403 Mon Sep 17 00:00:00 2001 From: Jack Chen Date: Mon, 8 Mar 2021 22:49:08 -0700 Subject: [PATCH 19/24] revised 3/9/2021 v4 --- .../utils/horizontal_interpolate.F90 | 20 +++++++------- src/chemistry/utils/tracer_data.F90 | 26 +++++++++---------- src/physics/cam/ssatcontrail.F90 | 4 +-- 3 files changed, 24 insertions(+), 26 deletions(-) diff --git a/src/chemistry/utils/horizontal_interpolate.F90 b/src/chemistry/utils/horizontal_interpolate.F90 index de0602213a..6a7ae7e4f8 100644 --- a/src/chemistry/utils/horizontal_interpolate.F90 +++ b/src/chemistry/utils/horizontal_interpolate.F90 @@ -111,7 +111,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight ! 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 ! |-------------------| @@ -122,7 +122,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight else weight_x(i2,i1) = (x1_east-x1_west)/(x2_east-x2_west) endif - elseif ( (x1_west.ge.x2_west).and.(x1_west.lt.x2_east) ) then + elseif ( (x1_west>=x2_west).and.(x1_westx2_west).and.(x1_east.le.x2_east) ) then + elseif ( (x1_east>x2_west).and.(x1_east<=x2_east) ) then ! case 3: ! x1_west x1_east ! |--------------------------------| @@ -144,7 +144,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight else weight_x(i2,i1) = (x1_east-x2_west)/(x2_east-x2_west) endif - elseif ( (x1_east.gt.x2_east).and.(x1_west.lt.x2_west) ) then + elseif ( (x1_east>x2_east).and.(x1_westslon2(im2+1)) then ! case 1: ! slon1(im1) slon1(im1+1) <--- end point ! |-------------------------| @@ -171,7 +171,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight 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 ! |-------------------| @@ -213,7 +213,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight else weight_y(j2,j1) = gw1(j1)/gw2(j2) endif - elseif ( (y1_south.ge.y2_south).and.(y1_south.lt.y2_north) ) then + elseif ( (y1_south>=y2_south).and.(y1_southy2_south).and.(y1_north<=y2_north) ) then ! case 3: ! y1_south y1_north ! |--------------------------------| @@ -235,7 +235,7 @@ subroutine xy_interp_init(im1,jm1,lon0,lat0,im2,jm2,weight_x,weight_y,use_flight else weight_y(j2,j1) = (y1_north-y2_south)/(y1_north-y1_south)*gw1(j1)/gw2(j2) endif - elseif ( (y1_north.gt.y2_north).and.(y1_south.lt.y2_south) ) then + elseif ( (y1_north>y2_north).and.(y1_south0.0_r8 ) then file%count_x(i2) = file%count_x(i2) + 1 file%index_x(i2,file%count_x(i2)) = i1 endif @@ -718,7 +718,7 @@ 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 @@ -1201,7 +1201,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 @@ -2389,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 @@ -2416,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 @@ -2475,7 +2475,7 @@ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, end if end do - if( tl.gt.src_x(nsrc+1)) sil = nsrc + if( tl>src_x(nsrc+1)) sil = nsrc y = 0.0_r8 bot = min(tl,src_x(nsrc+1)) @@ -2587,16 +2587,16 @@ subroutine vert_interp( ncol, levsiz, pin, pmid, datain, dataout ) ! do kk=kkstart,levsiz-1 do i=1,ncol - if (pin(i,kk).lt.pmid(i,k) .and. pmid(i,k).le.pin(i,kk+1)) then + if (pin(i,kk) 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/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index d2816f94b0..ed1e1695c0 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -7,11 +7,9 @@ module ssatcontrail 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_set_field, pbuf_old_tim_idx + 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 geopotential, only: geopotential_dse use phys_grid, only: get_wght_all_p - use time_manager, only: get_curr_date use wv_saturation, only: qsat_water, qsat_ice use aircraft_emit, only: get_aircraft From 49d07ae93c5a7052b835f12fabaa8743461813ad Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Tue, 9 Mar 2021 09:48:04 -0700 Subject: [PATCH 20/24] replace problematic tab --- src/chemistry/utils/tracer_data.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 6557c5b581..7615ef95c3 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -624,7 +624,7 @@ 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), stat=astat) + 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') From f8aa85ae59ce9db68cbb91bbc57554423f7f0238 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Wed, 17 Mar 2021 17:10:11 -0600 Subject: [PATCH 21/24] Fix some bugs and add cam regression test --- cime_config/testdefs/testlist_cam.xml | 9 ++++ .../cam/outfrq9s_contrail/shell_commands | 2 + .../cam/outfrq9s_contrail/user_nl_cam | 8 +++ .../cam/outfrq9s_contrail/user_nl_clm | 27 ++++++++++ src/chemistry/utils/tracer_data.F90 | 1 - src/physics/cam/ssatcontrail.F90 | 52 +++++++++---------- 6 files changed, 71 insertions(+), 28 deletions(-) create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/shell_commands create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_cam create mode 100644 cime_config/testdefs/testmods_dirs/cam/outfrq9s_contrail/user_nl_clm diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml index 66cf22b197..0ddf3034e3 100644 --- a/cime_config/testdefs/testlist_cam.xml +++ b/cime_config/testdefs/testlist_cam.xml @@ -1096,6 +1096,15 @@ + + + + + + + + + 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/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 7615ef95c3..6f4169bafc 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -1965,7 +1965,6 @@ subroutine interpolate_trcdata( state, flds, file, pbuf2d ) call vert_interp_mixrat(ncol,file%nlev,pver,state(c)%pint, & datain, data_out(:,:), & file%p0,ps,file%hyai,file%hybi,file%dist) - endif else call vert_interp(ncol, file%nlev, pin, state(c)%pmid, datain, data_out(:,:) ) endif diff --git a/src/physics/cam/ssatcontrail.F90 b/src/physics/cam/ssatcontrail.F90 index ed1e1695c0..17c620ae24 100644 --- a/src/physics/cam/ssatcontrail.F90 +++ b/src/physics/cam/ssatcontrail.F90 @@ -8,7 +8,7 @@ module ssatcontrail 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 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 @@ -23,7 +23,7 @@ module ssatcontrail 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) @@ -36,7 +36,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) !------------------------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 :: i,k integer :: lchnk,ncol real(r8) :: contrail(pcols,pver), pcontrail(pcols,pver) real(r8), pointer, dimension(:,:) :: cld ! cloud fraction @@ -49,12 +49,12 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) logical :: has_aircraft_H2O logical :: has_aircraft_distance real(r8) :: hkl, hkk, tv - real(r8) :: particle_mass + 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) :: dcld(pcols,pver) real(r8) :: ac_Q, ac_Q1, ac_Q2 real(r8) :: RHcts(pcols,pver) logical :: lq(pcnst) @@ -102,7 +102,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) 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 @@ -113,7 +113,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) lchnk = state1%lchnk ncol = state1%ncol - + call get_wght_all_p(lchnk, ncol, wght) itim = pbuf_old_tim_idx() @@ -125,9 +125,8 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) 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 + Ma = mwdry Mh2o = mwh2o ! contrail paramter (G = CF*p/epi) @@ -137,21 +136,20 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) 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(i,k) = 0._r8 - dcld(i,k) = 0._r8 + + 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 + + ! compute saturation pressure call qsat_water(T_contr, p, eslTc, qslTc) call qsat_water(state1%t(i,k), p, eslT, qslT) @@ -159,9 +157,9 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) ! 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) + + 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) @@ -170,7 +168,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) ! 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 +! 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)) @@ -181,7 +179,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) 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) + 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) @@ -192,10 +190,10 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) ! 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, + +! 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 @@ -203,7 +201,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) ! 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 @@ -221,7 +219,7 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) 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 + 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 @@ -230,5 +228,5 @@ subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc) enddo end subroutine ssatcontrail_d0 - + end module ssatcontrail From 3cc0e7fb50c6a7fe67f5108eeacc24b9378ea40b Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Fri, 19 Mar 2021 11:44:31 -0600 Subject: [PATCH 22/24] Remove extra test --- cime_config/testdefs/testlist_cam.xml | 1 - 1 file changed, 1 deletion(-) diff --git a/cime_config/testdefs/testlist_cam.xml b/cime_config/testdefs/testlist_cam.xml index 0ddf3034e3..7b91d76214 100644 --- a/cime_config/testdefs/testlist_cam.xml +++ b/cime_config/testdefs/testlist_cam.xml @@ -1099,7 +1099,6 @@ - From b1b1a31cabe65b10fc4c5dccad77a62c5f65efcf Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Mon, 22 Mar 2021 17:09:28 -0600 Subject: [PATCH 23/24] Fix MPI bug --- src/chemistry/utils/tracer_data.F90 | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index 6f4169bafc..957af638f1 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -114,7 +114,7 @@ module tracer_data 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 = .false. + logical :: dist real(r8) :: p0 type(var_desc_t) :: ps_id @@ -222,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' ) @@ -751,30 +752,30 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & endif endif - call mpi_bcast(file%weight_x, 1, mpi_real8 , mstrid, mpicom,ierr) + 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, 1, mpi_real8 , mstrid, mpicom,ierr) + 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, 1, mpi_integer , mstrid, mpicom,ierr) + 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, 1, mpi_integer , mstrid, mpicom,ierr) + 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, 1, mpi_integer , mstrid, mpicom,ierr) + 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, 1, mpi_integer , mstrid, mpicom,ierr) + 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, 1, mpi_real8 , mstrid, mpicom,ierr) + 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, 1, mpi_real8 , mstrid, mpicom,ierr) + 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, 1, mpi_integer , mstrid, mpicom,ierr) + 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, 1, mpi_integer , mstrid, mpicom,ierr) + 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, 1, mpi_integer , mstrid, mpicom,ierr) + 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, 1, mpi_integer , mstrid, mpicom,ierr) + 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 From 3340aba5fda549c33c6fbabdbb1846cd78ef5ef1 Mon Sep 17 00:00:00 2001 From: Cheryl Craig Date: Tue, 23 Mar 2021 13:48:40 -0600 Subject: [PATCH 24/24] Update ChangeLog --- doc/ChangeLog | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) 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