!CPP directives to control ic/bc conditions... !(The directive in module_mosaic_initmixrats also needs to be set.) ! CASENAME = 0 Uses Texas August 2004 case values and profiles ! 1 Uses same concentrations as TX, but uses different ! profiles depending on the species. (NEAQS2004 case) !#define CASENAME 0 MODULE module_input_smoke_data USE module_io_domain USE module_domain ! USE module_data_sorgam, ONLY : conmin, rgasuniv, epsilc, grav USE module_get_file_names, ONLY : eligible_file_name, number_of_eligible_files, unix_ls IMPLICIT NONE ! REAL :: grav = 9.8104 REAL, PARAMETER :: conmin=1.E-16, epsilc=1.E-16 ! Variables for adaptive time steps... TYPE(WRFU_Time), DIMENSION(max_domains) :: last_chem_time !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Initial atmospheric chemistry profile data INTEGER :: k_loop ! Used for loop index INTEGER :: lo ! number of chemicals in inital profile INTEGER :: logg ! number of final chemical species (nch-1) INTEGER :: kx ! number of vertical levels in temp profile INTEGER :: kxm1 PARAMETER( kx=16, kxm1=kx-1, logg=350, lo=34) ! DL (6/2/2013) Changed value of logg from 200 to 350 for additional gas species INTEGER, DIMENSION(logg) :: iref REAL, DIMENSION(logg) :: fracref REAL, DIMENSION(kx) :: dens REAL, DIMENSION(kx+1) :: zfa REAL, DIMENSION(kx+1) :: zfa_bdy REAL, DIMENSION(lo,kx) :: xl ! REAL :: so4vaptoaer ! DATA so4vaptoaer/.999/ ! CHARACTER (LEN=20), DIMENSION(logg) :: ggnam !BSINGH(12/04/13): changed length(LEN) from 4 to 20 ! type(lbc_concentration), allocatable :: fixed_lbc(:) CONTAINS SUBROUTINE vinterp_chem(nx1, nx2, ny1, ny2, nz1, nz_in, nz_out, nch, z_in, z_out, & data_in, data_out, extrapolate) ! Interpolates columns of chemistry data from one set of height surfaces to ! another INTEGER, INTENT(IN) :: nx1, nx2 INTEGER, INTENT(IN) :: ny1, ny2 INTEGER, INTENT(IN) :: nz1 INTEGER, INTENT(IN) :: nz_in INTEGER, INTENT(IN) :: nz_out INTEGER, INTENT(IN) :: nch REAL, INTENT(IN) :: z_in (nx1:nx2,nz1:nz_in ,ny1:ny2) REAL, INTENT(IN) :: z_out(nx1:nx2,nz1:nz_out,ny1:ny2) REAL, INTENT(IN) :: data_in (nx1:nx2,nz1:nz_in ,ny1:ny2,nch) REAL, INTENT(OUT) :: data_out(nx1:nx2,nz1:nz_out,ny1:ny2,nch) LOGICAL, INTENT(IN) :: extrapolate INTEGER :: i,j,l INTEGER :: k,kk REAL :: desired_z REAL :: dvaldz REAL :: wgt0 ! Loop over the number of chemical species chem_loop: DO l = 2, nch data_out(:,:,:,l) = -99999.9 DO j = ny1, ny2 DO i = nx1, nx2 output_loop: DO k = nz1, nz_out desired_z = z_out(i,k,j) IF (desired_z .LT. z_in(i,1,j)) THEN IF ((desired_z - z_in(i,1,j)).LT. 0.0001) THEN data_out(i,k,j,l) = data_in(i,1,j,l) ELSE IF (extrapolate) THEN ! Extrapolate downward because desired height level is below ! the lowest level in our input data. Extrapolate using simple ! 1st derivative of value with respect to height for the bottom 2 ! input layers. ! Add a check to make sure we are not using the gradient of ! a very thin layer IF ( (z_in(i,1,j) - z_in(i,2,j)) .GT. 0.001) THEN dvaldz = (data_in(i,1,j,l) - data_in(i,2,j,l)) / & (z_in(i,1,j) - z_in(i,2,j) ) ELSE dvaldz = (data_in(i,1,j,l) - data_in(i,3,j,l)) / & (z_in(i,1,j) - z_in(i,3,j) ) ENDIF data_out(i,k,j,l) = MAX( data_in(i,1,j,l) + & dvaldz * (desired_z-z_in(i,1,j)), 0.) ELSE data_out(i,k,j,l) = data_in(i,1,j,l) ENDIF ENDIF ELSE IF (desired_z .GT. z_in(i,nz_in,j)) THEN IF ( (z_in(i,nz_in,j) - desired_z) .LT. 0.0001) THEN data_out(i,k,j,l) = data_in(i,nz_in,j,l) ELSE IF (extrapolate) THEN ! Extrapolate upward IF ( (z_in(i,nz_in-1,j)-z_in(i,nz_in,j)) .GT. 0.0005) THEN dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-1,j,l)) / & (z_in(i,nz_in,j) - z_in(i,nz_in-1,j)) ELSE dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-2,j,l)) / & (z_in(i,nz_in,j) - z_in(i,nz_in-2,j)) ENDIF data_out(i,k,j,l) = MAX( data_in(i,nz_in,j,l) + & dvaldz * (desired_z-z_in(i,nz_in,j)), 0.) ELSE data_out(i,k,j,l) = data_in(i,nz_in,j,l) ENDIF ENDIF ELSE ! We can trap between two levels and linearly interpolate input_loop: DO kk = 1, nz_in-1 IF (desired_z .EQ. z_in(i,kk,j) )THEN data_out(i,k,j,l) = data_in(i,kk,j,l) EXIT input_loop ELSE IF (desired_z .EQ. z_in(i,kk+1,j) )THEN data_out(i,k,j,l) = data_in(i,kk+1,j,l) EXIT input_loop ELSE IF ( (desired_z .GT. z_in(i,kk,j)) .AND. & (desired_z .LT. z_in(i,kk+1,j)) ) THEN wgt0 = (desired_z - z_in(i,kk+1,j)) / & (z_in(i,kk,j)-z_in(i,kk+1,j)) data_out(i,k,j,l) = MAX( wgt0*data_in(i,kk,j,l) + & (1.-wgt0)*data_in(i,kk+1,j,l), 0.) EXIT input_loop ENDIF ENDDO input_loop ENDIF ENDDO output_loop ENDDO ENDDO ENDDO chem_loop RETURN END SUBROUTINE vinterp_chem !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE input_chem_profile (si_grid) IMPLICIT NONE TYPE(domain) :: si_grid INTEGER :: i , j , k, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe INTEGER :: fid, ierr, numgas INTEGER :: debug_level REAL, ALLOCATABLE, DIMENSION(:,:,:) :: si_zsigf, si_zsig CHARACTER (LEN=80) :: inpname, message write(message,'(A)') 'Subroutine input_chem_profile: ' CALL wrf_message ( TRIM(message) ) ! And here is an instance of using the information in the NAMELIST. CALL nl_get_debug_level ( 1,debug_level ) CALL set_wrf_debug_level ( debug_level ) ! Get grid dimensions CALL get_ijk_from_grid ( si_grid , & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe ) IF ( si_grid%chem_opt == CHEM_TRACER ) THEN si_grid%chem(ims:ime,kms:kme,jms:jme,:) = 0.0 ENDIF RETURN END SUBROUTINE input_chem_profile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE make_chem_profile ( nx1, nx2, ny1, ny2, nz1, nz2, nch, numgas, & chem_opt, zgrid, chem ) IMPLICIT NONE INTEGER, INTENT(IN) :: nx1, ny1, nz1 INTEGER, INTENT(IN) :: nx2, ny2, nz2 INTEGER, INTENT(IN) :: nch, numgas, chem_opt !REAL, INTENT(IN), DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2) :: zgrid REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2) :: zgrid CHARACTER (LEN=80) :: message INTEGER :: i, j, k, l, is REAL, DIMENSION(nx1:nx2,nz1:kx ,ny1:ny2,lo+1):: chprof REAL, DIMENSION(nx1:nx2,nz1:kx ,ny1:ny2) :: zprof REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2,nch) :: chem REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2,lo ) :: stor REAL :: hc358 !wig, 2-May-2007 REAL :: olit !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Check the number of species if( nch .NE. num_chem) then message = ' Input_chem_profile: wrong number of chemical species' ! CALL WRF_ERROR_FATAL ( message ) endif ! Vertically flip the chemistry data as it is given top down and ! heights are bottom up. Fill temp 3D chemical and profile array, ! keep chem slot 1 open as vinterp_chem assumes there is no data. DO j=ny1,ny2 DO k= 1,kx DO i=nx1,nx2 chprof(i,k,j,2:lo+1) = xl(1:lo,kx-k+1) zprof(i,k,j) = 0.5*(zfa(k)+zfa(k+1)) ENDDO ENDDO ENDDO ! ! return xl to previous value for next time... ! 34 chemicals (lo), 16 vertical levels (kx) ! DO i=lo-6,lo ! xl(i,1:kx)=xl(i,1:kx)*dens(1:kx) ! ENDDO ! Change number concentrations to mixing ratios for short-lived NALROM species do k=1,kx chprof(:,k,:,lo-5:lo+1) = chprof(:,k,:,lo-5:lo+1)/dens(k) end do ! Interpolate temp 3D chemical and profile array to WRF grid call vinterp_chem(nx1, nx2, ny1, ny2, nz1, kx, nz2, lo, zprof, zgrid, & chprof, chem, .false.) ! place interpolated data into temp storage array stor(nx1:nx2,nz1:nz2,ny1:ny2,1:lo) = chem(nx1:nx2,nz1:nz2,ny1:ny2,2:lo+1) ! Here is where the chemistry profile is constructed !chem(:,:,:,1) = stor(:,:,:,1) * 0. chem(nx1:nx2,nz1:nz2,ny1:ny2,1) = -999. ! DO l=2,nch DO l=2,numgas is=iref(l-1) DO j=ny1,ny2 DO k=nz1,nz2 DO i=nx1,nx2 chem(i,k,j,l)=fracref(l-1)*stor(i,k,j,is)*1.E6 ENDDO ENDDO ENDDO ENDDO ! ! For CBMZ, we need to construct PAR based on a combination of other ! species. This cannot be done with the looping construct above so ! we have to treat it separately. (wig, 2-May-2007) ! SELECT CASE(chem_opt) ! CASE (CBMZ,CBMZ_BB,CBMZ_BB_KPP, CBMZ_MOSAIC_KPP, & ! CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN, & ! CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ, & ! CBMZSORG,CBMZSORG_AQ, CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, & ! CBMZ_MOSAIC_DMS_4BIN_AQ, CBMZ_MOSAIC_DMS_8BIN_AQ, & ! CBMZ_CAM_MAM3_NOAQ, CBMZ_CAM_MAM3_AQ, CBMZ_CAM_MAM7_NOAQ, CBMZ_CAM_MAM7_AQ) ! do j = ny1,ny2 ! do k = nz1,nz2 ! do i = nx1,nx2 !Construct the sum of the profiles for hc3, hc5, & hc8 ! hc358 = ( 2.9*fracref(numgas+1)*stor(i,k,j,iref(numgas+1)) & ! +4.8*fracref(numgas+2)*stor(i,k,j,iref(numgas+2)) & ! +7.9*fracref(numgas+3)*stor(i,k,j,iref(numgas+3)) & ! )*1.E6 ! chem(i,k,j,p_par) = & ! 0.4*chem(i,k,j,p_ald) + hc358 & ! + 0.9*chem(i,k,j,p_ket) + 2.8*chem(i,k,j,p_oli) & ! + 1.8*chem(i,k,j,p_olt) + 1.0*chem(i,k,j,p_ora2) ! end do ! end do ! end do CASE (CBM4_KPP) do j = ny1,ny2 do k = nz1,nz2 do i = nx1,nx2 !Construct the sum of the profiles for hc3, hc5, & hc8 hc358 = ( 2.9*fracref(numgas+1)*stor(i,k,j,iref(numgas+1)) & +4.8*fracref(numgas+2)*stor(i,k,j,iref(numgas+2)) & +7.9*fracref(numgas+3)*stor(i,k,j,iref(numgas+3)) & )*1.E6 olit = ( 0.9*fracref(numgas+4)*stor(i,k,j,iref(numgas+4)) & +2.8*fracref(numgas+5)*stor(i,k,j,iref(numgas+5)) & +1.8*fracref(numgas+6)*stor(i,k,j,iref(numgas+6)) & +1.0*fracref(numgas+7)*stor(i,k,j,iref(numgas+7)) & )*1.E6 chem(i,k,j,p_par) = 0.4*chem(i,k,j,p_ald2) + hc358 + olit end do end do end do CASE (CB05_SORG_AQ_KPP) do j = ny1,ny2 do k = nz1,nz2 do i = nx1,nx2 !Construct the sum of the profiles for hc3, hc5, & hc8 hc358 = ( 2.9*fracref(numgas+1)*stor(i,k,j,iref(numgas+1)) & +4.8*fracref(numgas+2)*stor(i,k,j,iref(numgas+2)) & +7.9*fracref(numgas+3)*stor(i,k,j,iref(numgas+3)) & )*1.E6 chem(i,k,j,p_par) = & 0.4*chem(i,k,j,p_ald2) + hc358 & +0.4*chem(i,k,j,p_aldx) + 2.8*chem(i,k,j,p_ole) & + 1.8*chem(i,k,j,p_iole) + 1.0*chem(i,k,j,p_aacd) end do end do end do CASE (CB05_SORG_VBS_AQ_KPP) do j = ny1,ny2 do k = nz1,nz2 do i = nx1,nx2 !Construct the sum of the profiles for hc3, hc5, & hc8 hc358 = ( 2.9*fracref(numgas+1)*stor(i,k,j,iref(numgas+1)) & +4.8*fracref(numgas+2)*stor(i,k,j,iref(numgas+2)) & +7.9*fracref(numgas+3)*stor(i,k,j,iref(numgas+3)) & )*1.E6 chem(i,k,j,p_par) = & 0.4*chem(i,k,j,p_ald2) + hc358 & +0.4*chem(i,k,j,p_aldx) + 2.8*chem(i,k,j,p_ole) & + 1.8*chem(i,k,j,p_iole) + 1.0*chem(i,k,j,p_aacd) end do end do end do END SELECT RETURN END SUBROUTINE make_chem_profile !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE bdy_chem_value_tracer ( chem, nch ) ! This subroutine is called to set the boundary values of chemistry ! species when chem_opt==CHEM_TRACER. Typically, the boundary values ! here should be set to match those in input_chem_profile so that the ! interior and boundary values are the same. ! William.Gustafson@pnl.gov; 16-Jun-2005 IMPLICIT NONE REAL, intent(OUT) :: chem INTEGER, intent(IN) :: nch ! index number of chemical species !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if( nch .ne. p_co ) then chem = 0.0001 else if( nch == p_co ) then chem = 0.08 else chem = conmin end if if( nch .eq. p_tracer_1 ) then chem = 0.08 endif END SUBROUTINE bdy_chem_value_tracer !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE bdy_chem_value ( chem, z, nch, numgas ) IMPLICIT NONE REAL, intent(OUT) :: chem REAL, intent(IN) :: z ! 3D height array INTEGER, intent(IN) :: nch ! index number of chemical species INTEGER, intent(IN) :: numgas ! index number of last gas species INTEGER :: i, k, irefcur REAL, DIMENSION(kx):: cprof ! chemical profile, diff. index order REAL, DIMENSION(1:kx):: zprof REAL :: stor REAL :: wgt0 CHARACTER (LEN=80) :: message !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Check the number of species ! if((nch-1).gt.logg)return ! for radmkpp there is co2 and ch4 in the variable list ! if(p_co2.gt.1)then if (nch.eq.p_co2)then chem=370. return endif if (nch.eq.p_ch4)then chem=1.7 return endif endif ! if( nch .GT. numgas) then ! message = ' Input_chem_profile: wrong number of chemical species' ! return ! CALL WRF_ERROR_FATAL ( message ) ! endif ! Vertically flip the chemistry data as it is given top down ! and heights in zfa are bottom up ! Fill 1D chemical profile array cprof ! Convert species 28-34 (lo-6:lo) from (molecules/cm3) to (mol/mol) irefcur = iref(nch-1) DO k = 1,kx zprof(k) = 0.5*(zfa_bdy(k)+zfa_bdy(k+1)) if (irefcur .lt. lo-6) then cprof(k) = xl(irefcur,kx+1-k) else cprof(k) = xl(irefcur,kx+1-k)/dens(kx+1-k) end if ENDDO ! Interpolate temp 3D chemical profile array to WRF field IF (z .LT. zprof(1)) THEN stor = cprof(1) ELSE IF (z .GT. zprof(kx)) THEN stor = cprof(kx) ELSE ! We can trap between two levels and linearly interpolate input_loop: DO k = 1, kx-1 IF (z .EQ. zprof(k) )THEN stor = cprof(k) EXIT input_loop ELSE IF (z .EQ. zprof(k+1) )THEN stor = cprof(k+1) EXIT input_loop ELSE IF ( (z .GT. zprof(k)) .AND. & (z .LT. zprof(k+1)) ) THEN wgt0 = (z - zprof(k+1)) / & (zprof(k) - zprof(k+1)) stor = MAX( wgt0 *cprof(k ) + & (1.-wgt0)*cprof(k+1), 0.) EXIT input_loop ENDIF ENDDO input_loop ENDIF ! Here is where the chemistry value is constructed chem = fracref(nch-1)*stor*1.E6 ! special code for sulfate/h2so4 ! if(nch.eq.p_sulf.and.p_nu0.gt.1)then ! chem=chem*(1.-so4vaptoaer) ! endif RETURN END SUBROUTINE bdy_chem_value !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE flow_dep_bdy_smoke (chem, & chem_bxs,chem_btxs, & chem_bxe,chem_btxe, & chem_bys,chem_btys, & chem_bye,chem_btye, & dtbdy, & !have_bcs_chem, & u, v, config_flags, alt, & ! t,pb,p, & ic, & ids,ide, jds,jde, kds,kde, & ! domain dims ims,ime, jms,jme, kms,kme, & ! memory dims ips,ipe, jps,jpe, kps,kpe, & ! patch dims its,ite, jts,jte, kts,kte ) ! This subroutine is called regardless of have_bcs_chem value in the namelist. ! This subroutine sets zero gradient conditions for outflow and a set profile value ! for inflow in the boundary specified region. Note that field must be unstaggered. ! The velocities, u and v, will only be used to check their sign (coupled vels OK) ! spec_zone is the width of the outer specified b.c.s that are set here. ! (JD August 2000) ! USE module_cam_mam_initmixrats, only: bdy_chem_value_cam_mam IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme INTEGER, INTENT(IN ) :: ips,ipe, jps,jpe, kps,kpe INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte INTEGER, INTENT(IN ) :: ic REAL, INTENT(IN ) :: dtbdy TYPE( grid_config_rec_type ) :: config_flags REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: chem REAL, DIMENSION( jms:jme , kds:kde , config_flags%spec_bdy_width), INTENT(IN ) :: chem_bxs, chem_bxe, chem_btxs, chem_btxe REAL, DIMENSION( ims:ime , kds:kde , config_flags%spec_bdy_width), INTENT(IN ) :: chem_bys, chem_bye, chem_btys, chem_btye ! REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: z REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alt REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: u REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: v ! real, INTENT (IN) :: rcp,t0,p1000mb INTEGER :: i, j, k, numgas INTEGER :: ibs, ibe, jbs, jbe, itf, jtf, ktf INTEGER :: i_inner, j_inner INTEGER :: b_dist integer :: itestbc, i_bdy_method, spec_zone, spec_bdy_width ! real tempfac,convfac real :: chem_bv_def logical :: have_bcs_chem INTEGER, SAVE :: icall chem_bv_def = conmin numgas = 0 !get_last_gas(config_flags%chem_opt) itestbc=0 ibs = ids ibe = ide-1 itf = min(ite,ide-1) jbs = jds jbe = jde-1 jtf = min(jte,jde-1) ktf = kde-1 have_bcs_chem = config_flags%have_bcs_chem spec_bdy_width= config_flags%spec_bdy_width spec_zone = config_flags%spec_zone i_bdy_method = 0 if (have_bcs_chem) i_bdy_method =6 if (ic .lt. param_first_scalar) i_bdy_method = 0 IF (jts - jbs .lt. spec_zone) THEN ! Y-start boundary, first boundary in south-north direction, j index is for SN direction, i index is for WE direction DO j = jts, min(jtf,jbs+spec_zone-1) b_dist = j - jbs DO k = kts, ktf DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist) i_inner = max(i,ibs+spec_zone) i_inner = min(i_inner,ibe-spec_zone) IF(v(i,k,j) .lt. 0.)THEN chem(i,k,j) = chem(i_inner,k,jbs+spec_zone) ! air parcel leaving the domain ELSE ! RAR: if (i_bdy_method .eq. 6) then ! have_bcs_chem=.true. CALL bdy_chem_value_rapsmoke( chem(i,k,j),chem_bys(i,k,1),chem_btys(i,k,1),alt(i,k,j),dtbdy ) ! air parcel entering the domain else chem(i,k,j) = chem_bv_def endif ENDIF ENDDO ENDDO ENDDO ENDIF IF (jbe - jtf .lt. spec_zone) THEN ! Y-end boundary DO j = max(jts,jbe-spec_zone+1), jtf b_dist = jbe - j DO k = kts, ktf DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist) i_inner = max(i,ibs+spec_zone) i_inner = min(i_inner,ibe-spec_zone) IF(v(i,k,j+1) .gt. 0.)THEN chem(i,k,j) = chem(i_inner,k,jbe-spec_zone) ELSE ! if (i_bdy_method .eq. 1) then ! CALL bdy_chem_value ( & ! chem(i,k,j), z(i,k,j), ic, numgas ) ! else if (i_bdy_method .eq. 9) then ! CALL bdy_chem_value_racm ( & ! chem(i,k,j), z(i,k,j), ic, numgas,p_co2 ) ! else if (i_bdy_method .eq. 2) then ! tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp ! convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac ! CALL bdy_chem_value_sorgam ( & ! chem(i,k,j), z(i,k,j), ic, config_flags, & ! alt(i,k,j),convfac,g) ! else if (i_bdy_method .eq. 3) then ! CALL bdy_chem_value_cbmz ( & ! chem(i,k,j), z(i,k,j), ic, numgas ) ! else if (i_bdy_method .eq. 4) then ! CALL bdy_chem_value_mosaic ( & ! chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags ) ! else if (i_bdy_method .eq. 5) then ! CALL bdy_chem_value_tracer ( chem(i,k,j), ic ) ! RAR: if (i_bdy_method .eq. 6) then ! have_bcs_chem=.true. CALL bdy_chem_value_rapsmoke ( chem(i,k,j),chem_bye(i,k,1),chem_btye(i,k,1),alt(i,k,j),dtbdy) ! else if (i_bdy_method .eq. 7) then ! CALL bdy_chem_value_gocart ( chem(i,k,j), ic ) ! else if (i_bdy_method .eq. 8) then ! chem(i,k,j) = 0. ! else if (i_bdy_method .eq. 16) then ! CALL bdy_chem_value_ghg ( chem(i,k,j), ic ) ! For GHGs ! else if (i_bdy_method .eq. 15) then ! CALL bdy_chem_value_cb05 ( & ! 2, chem(i,k,j), k, ic, config_flags, numgas ) ! else if (i_bdy_method .eq. 17) then ! CALL bdy_chem_value_cb05_vbs ( & ! 2, chem(i,k,j), k, ic, config_flags, numgas ) ! else if (i_bdy_method .eq. 501) then ! tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp ! convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac ! CALL bdy_chem_value_cam_mam( & ! chem(i,k,j), z(i,k,j), ic, config_flags, alt(i,k,j), convfac, g ) else chem(i,k,j) = chem_bv_def endif ENDIF ENDDO ENDDO ENDDO ENDIF IF (its - ibs .lt. spec_zone) THEN ! X-start boundary DO i = its, min(itf,ibs+spec_zone-1) b_dist = i - ibs DO k = kts, ktf DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1) j_inner = max(j,jbs+spec_zone) j_inner = min(j_inner,jbe-spec_zone) IF(u(i,k,j) .lt. 0.)THEN chem(i,k,j) = chem(ibs+spec_zone,k,j_inner) ELSE ! if (i_bdy_method .eq. 1) then ! CALL bdy_chem_value ( & ! chem(i,k,j), z(i,k,j), ic, numgas ) ! else if (i_bdy_method .eq. 9) then ! CALL bdy_chem_value_racm ( & ! chem(i,k,j), z(i,k,j), ic, numgas,p_co2 ) ! else if (i_bdy_method .eq. 2) then ! tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp ! convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac ! CALL bdy_chem_value_sorgam ( & ! chem(i,k,j), z(i,k,j), ic, config_flags, & ! alt(i,k,j),convfac,g) ! else if (i_bdy_method .eq. 3) then ! CALL bdy_chem_value_cbmz ( & ! chem(i,k,j), z(i,k,j), ic, numgas ) ! else if (i_bdy_method .eq. 4) then ! CALL bdy_chem_value_mosaic ( & ! chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags ) ! else if (i_bdy_method .eq. 5) then ! CALL bdy_chem_value_tracer ( chem(i,k,j), ic ) !RAR: west boundary in WE direction if (i_bdy_method .eq. 6) then ! have_bcs_chem=.true. CALL bdy_chem_value_rapsmoke (chem(i,k,j),chem_bxs(j,k,1),chem_btxs(j,k,1),alt(i,k,j),dtbdy) ! else if (i_bdy_method .eq. 7) then ! CALL bdy_chem_value_gocart ( chem(i,k,j), ic ) ! else if (i_bdy_method .eq. 8) then ! chem(i,k,j) = 0. ! else if (i_bdy_method .eq. 16) then ! CALL bdy_chem_value_ghg ( chem(i,k,j), ic ) ! For GHGs ! else if (i_bdy_method .eq. 15) then ! CALL bdy_chem_value_cb05 ( & ! 3, chem(i,k,j), k, ic, config_flags, numgas ) ! else if (i_bdy_method .eq. 17) then ! CALL bdy_chem_value_cb05_vbs ( & ! 3, chem(i,k,j), k, ic, config_flags, numgas ) ! else if (i_bdy_method .eq. 501) then ! tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp ! convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac ! CALL bdy_chem_value_cam_mam( & ! chem(i,k,j), z(i,k,j), ic, config_flags, alt(i,k,j), convfac, g ) else chem(i,k,j) = chem_bv_def endif ENDIF ENDDO ENDDO ENDDO ENDIF IF (ibe - itf .lt. spec_zone) THEN ! X-end boundary DO i = max(its,ibe-spec_zone+1), itf b_dist = ibe - i DO k = kts, ktf DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1) j_inner = max(j,jbs+spec_zone) j_inner = min(j_inner,jbe-spec_zone) IF(u(i+1,k,j) .gt. 0.)THEN chem(i,k,j) = chem(ibe-spec_zone,k,j_inner) ELSE ! if (i_bdy_method .eq. 1) then ! CALL bdy_chem_value ( & ! chem(i,k,j), z(i,k,j), ic, numgas ) ! else if (i_bdy_method .eq. 9) then ! CALL bdy_chem_value_racm ( & ! chem(i,k,j), z(i,k,j), ic, numgas,p_co2 ) ! else if (i_bdy_method .eq. 2) then ! tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp ! convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac ! CALL bdy_chem_value_sorgam ( & ! chem(i,k,j), z(i,k,j), ic, config_flags, & ! alt(i,k,j),convfac,g) ! else if (i_bdy_method .eq. 3) then ! CALL bdy_chem_value_cbmz ( & ! chem(i,k,j), z(i,k,j), ic, numgas ) ! else if (i_bdy_method .eq. 4) then ! CALL bdy_chem_value_mosaic ( & ! chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags ) ! else if (i_bdy_method .eq. 5) then ! CALL bdy_chem_value_tracer ( chem(i,k,j), ic ) ! RAR: if (i_bdy_method .eq. 6) then ! have_bcs_chem=.true. CALL bdy_chem_value_rapsmoke ( chem(i,k,j),chem_bxe(j,k,1),chem_btxe(j,k,1),alt(i,k,j),dtbdy) ! else if (i_bdy_method .eq. 7) then ! CALL bdy_chem_value_gocart ( chem(i,k,j), ic ) ! else if (i_bdy_method .eq. 8) then ! chem(i,k,j) = 0. ! else if (i_bdy_method .eq. 16) then ! CALL bdy_chem_value_ghg ( chem(i,k,j), ic ) ! For GHGs ! else if (i_bdy_method .eq. 15) then ! CALL bdy_chem_value_cb05 ( & ! 4, chem(i,k,j), k, ic, config_flags, numgas ) ! else if (i_bdy_method .eq. 17) then ! CALL bdy_chem_value_cb05_vbs ( & ! 4, chem(i,k,j), k, ic, config_flags, numgas ) ! else if (i_bdy_method .eq. 501) then ! tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp ! convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac ! CALL bdy_chem_value_cam_mam( & ! chem(i,k,j), z(i,k,j), ic, config_flags, alt(i,k,j), convfac, g ) else chem(i,k,j) = chem_bv_def endif ENDIF ENDDO ENDDO ENDDO ENDIF IF (config_flags%debug_chem .AND. icall<2000) THEN WRITE(*,*) 'flow_dep_bdy_smoke: dtbdy=dt_rk+dtbc ',dtbdy WRITE(*,*) 'flow_dep_bdy_smoke: its,ite,jts,jte,kts,kte ',its,ite,jts,jte,kts,kte WRITE(*,*) 'flow_dep_bdy_smoke: alt(its,kts,jts),alt(ite,kte,jte) ',alt(its,kts,jts),alt(ite,kte,jte) WRITE(*,*) 'flow_dep_bdy_smoke: chem(its,kts,jts),chem_bxe(jts,kts,1),chem_btxe(jts,kts,1) ',chem(its,kts,jts),chem_bxe(jts,kts,1),chem_btxe(jts,kts,1) icall=icall+1 END IF END SUBROUTINE flow_dep_bdy_smoke !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE bdy_chem_value_rapsmoke ( chem, chem_b, chem_bt, irho, dt) IMPLICIT NONE REAL, intent(OUT) :: chem REAL, intent(IN) :: chem_b REAL, intent(IN) :: chem_bt, irho ! inverse density REAL, intent(IN) :: dt CHARACTER (LEN=80) :: message ! This field is in ug/m3 in the met_em files generated from RAP output, therefore we divide it here by air den chem=max(epsilc,chem_b + chem_bt * dt) * irho RETURN END SUBROUTINE bdy_chem_value_rapsmoke !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE cv_mmdd_jday ( YY, MM, DD, JDAY) ! ! Subroutine to compute the julian day given the month and day ! ! INTEGER, INTENT(IN ) :: YY, MM, DD INTEGER, INTENT(OUT) :: JDAY INTEGER, DIMENSION(12) :: imon, imon_a INTEGER :: i DATA imon_a /0,31,59,90,120,151,181,212,243,273,304,334/ ! !..... Check for leap year. ! do i=1,12 imon(i) = imon_a(i) enddo if(YY .eq. (YY/4)*4) then do i=3,12 imon(i) = imon(i) + 1 enddo endif ! !..... Convert month, day to julian day. ! jday = imon(mm) + dd END SUBROUTINE cv_mmdd_jday !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! integer FUNCTION get_last_gas(chem_opt) implicit none integer, intent(in) :: chem_opt ! Determine the index of the last gas species, which depends ! upon the gas mechanism. select case (chem_opt) case (0) get_last_gas = 0 !!! TUCCELLA case (RADM2, RADM2_KPP, RADM2SORG, RADM2SORG_AQ, RADM2SORG_AQCHEM, RADM2SORG_KPP, & RACM_KPP, RACMPM_KPP, RACM_MIM_KPP, RACMSORG_AQ, RACMSORG_AQCHEM_KPP, & RACM_ESRLSORG_AQCHEM_KPP, RACM_ESRLSORG_KPP, RACMSORG_KPP, RACM_SOA_VBS_KPP,& RACM_SOA_VBS_AQCHEM_KPP,GOCARTRACM_KPP,GOCARTRADM2) get_last_gas = p_ho2 ! case (SAPRC99_KPP,SAPRC99_MOSAIC_4BIN_VBS2_KPP, & ! SAPRC99_MOSAIC_8BIN_VBS2_AQ_KPP,SAPRC99_MOSAIC_8BIN_VBS2_KPP )!BSINGH(12/13/2013): Added SAPRC 8 bin AQ case and non-aq on 04/03/2014 ! get_last_gas = p_ch4 ! case (CBMZ,CBMZ_MOSAIC_DMS_4BIN,CBMZ_MOSAIC_DMS_8BIN,CBMZ_MOSAIC_DMS_4BIN_AQ,CBMZ_MOSAIC_DMS_8BIN_AQ) ! get_last_gas = p_mtf ! case (CBMZ_BB,CBMZ_BB_KPP, CBMZ_MOSAIC_KPP, CBMZ_MOSAIC_4BIN, & ! CBMZ_MOSAIC_8BIN,CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ,CBMZSORG,CBMZSORG_AQ) ! get_last_gas = p_isopo2 case (CHEM_TRACER) get_last_gas = p_co case (CHEM_TRACE2) get_last_gas = p_tracer_1 case (GOCART_SIMPLE) get_last_gas = p_msa case (CBM4_KPP) get_last_gas = p_ho2 case (CB05_SORG_AQ_KPP, CB05_SORG_VBS_AQ_KPP) get_last_gas = p_nh3 case (CHEM_VASH) get_last_gas = 0 case (CHEM_VOLC) get_last_gas = p_sulf case (CHEM_VOLC_4BIN) get_last_gas = 0 case (DUST) get_last_gas = 0 case (MOZART_KPP) get_last_gas = p_meko2 ! case (CRIMECH_KPP, CRI_MOSAIC_8BIN_AQ_KPP, CRI_MOSAIC_4BIN_AQ_KPP) ! GET_LAST_GAS = p_ic3h7no3 ! case (MOZCART_KPP) ! get_last_gas = p_meko2 ! case (MOZART_MOSAIC_4BIN_KPP) ! get_last_gas = p_meko2 ! case (MOZART_MOSAIC_4BIN_AQ_KPP) ! get_last_gas = p_meko2 case (CO2_TRACER,GHG_TRACER) ! No gas chemistry or deposition for GHGs get_last_gas = 0 case (CHEM_SMOKE) get_last_gas = 0 case ( CBMZ_CAM_MAM3_NOAQ, CBMZ_CAM_MAM3_AQ, CBMZ_CAM_MAM7_NOAQ, CBMZ_CAM_MAM7_AQ ) get_last_gas = p_soag case default call wrf_error_fatal("get_last_gas: could not decipher chem_opt value") end select END FUNCTION get_last_gas !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE sorgam_set_aer_bc_pnnl( chem, z, nch, config_flags ) ! USE module_data_sorgam, ONLY : dginia, dginin, dginic, esn36, esc36, esa36, seasfac, no3fac, nh4fac, so4fac, soilfac, anthfac, orgfac implicit none INTEGER,INTENT(IN ) :: nch real,intent(in ) :: z REAL,INTENT(INOUT ) :: chem TYPE (grid_config_rec_type) , INTENT (in) :: config_flags REAL :: mult, & m3acc, m3cor, m3nuc, & bv_so4ai, bv_so4aj, & bv_nh4ai, bv_nh4aj, & bv_no3ai, bv_no3aj, & bv_eci, bv_ecj, & bv_p25i, bv_p25j, & bv_orgpai,bv_orgpaj, & bv_antha, bv_seas, bv_soila ! ! Determine height multiplier... ! This should mimic the calculation in sorgam_init_aer_ic_pnnl, ! mosaic_init_wrf_mixrats_opt2, and bdy_chem_value_mosaic !!$! Height(m) Multiplier !!$! --------- ---------- !!$! <=2000 1.0 !!$! 2000=5000 0.25 !!$! !!$! which translates to: !!$! 2000 2000. & !!$ .and. z <= 3000. ) then !!$ mult = 1.0 - 0.0005*(z-2000.) !!$ elseif( z > 3000. & !!$ .and. z <= 5000. ) then !!$ mult = 0.5 - 1.25e-4*(z-3000.) !!$ else !!$ mult = 0.25 !!$ end if ! Updated aerosol profile multiplier 1-Apr-2005: ! Height(m) Multiplier ! --------- ---------- ! <=2000 1.0 ! 2000=5000 0.125 ! ! which translates to: ! 2000 2000. & ! .and. z <= 3000. ) then ! mult = 1.0 - 0.00075*(z-2000.) ! elseif( z > 3000. & ! .and. z <= 5000. ) then ! mult = 0.25 - 4.166666667e-5*(z-3000.) ! else ! mult = 0.125 ! end if ! if( z <= 500. ) then ! mult = 1.0 ! elseif( z > 500. & ! .and. z <= 1000. ) then ! mult = 1.0 - 0.001074*(z-500.) ! elseif( z > 1000. & ! .and. z <= 5000. ) then ! mult = 0.463 - 0.000111*(z-1000.) ! else ! mult = 0.019 ! end if ! These should match what is in sorgam_init_aer_ic_pnnl. ! Boundary values as of 2-Dec-2004: ! bv_so4aj = mult*2.375 ! bv_so4ai = mult*0.179 ! bv_nh4aj = mult*0.9604 ! bv_nh4ai = mult*0.0196 ! bv_no3aj = mult*0.0650 ! bv_no3ai = mult*0.0050 ! bv_ecj = mult*0.1630 ! bv_eci = mult*0.0120 ! bv_p25j = mult*0.6350 ! bv_p25i = mult*0.0490 ! bv_orgpaj = mult*0.9300 ! bv_orgpai = mult*0.0700 ! bv_antha = mult*2.2970 ! bv_seas = mult*0.2290 ! bv_soila = conmin !#if (CASENAME == 4) ! if( z <= 2000. ) then ! mult = 1.0 ! elseif( z > 2000. & ! .and. z <= 3000. ) then ! mult = 1.0 - 0.00075*(z-2000.) ! elseif( z > 3000. & ! .and. z <= 5000. ) then ! mult = 0.25 - 4.166666667e-5*(z-3000.) ! else ! mult = 0.125 ! end if ! bv_so4aj = mult*(0.0004810001+0.7271175)*0.97 ! bv_so4ai = mult*(0.0004810001+0.7271175)*0.03 ! bv_nh4aj = mult*0.2133708*0.97 ! bv_nh4ai = mult*0.2133708*0.03 ! bv_no3aj = mult*0.01399485*0.97 ! bv_no3ai = mult*0.01399485*0.03 ! bv_ecj = mult*0.04612048*0.97 ! bv_eci = mult*0.04612048*0.03 ! bv_p25j = mult*1.890001e-05*0.97 ! bv_p25i = mult*1.890001e-05*0.03 ! bv_antha = conmin ! bv_orgpaj = mult*0.5844942*0.97 ! bv_orgpai = mult*0.5844942*0.03 ! bv_seas = conmin ! bv_soila = conmin !#endif ! m3... calculations should match the very end of module_aerosols_sorgam.F !... i-mode (note that the 8 SOA species have bv=conmin) ! m3nuc = so4fac*bv_so4ai + nh4fac*bv_nh4ai + & ! no3fac*bv_no3ai + & ! orgfac*8.0*conmin + orgfac*bv_orgpai + & ! anthfac*bv_p25i + anthfac*bv_eci !... j-mode (note that the 8 SOA species have bv=conmin) ! m3acc = so4fac*bv_so4aj + nh4fac*bv_nh4aj + & ! no3fac*bv_no3aj + & ! orgfac*8.0*conmin + orgfac*bv_orgpaj + & ! anthfac*bv_p25j + anthfac*bv_ecj !...c-mode ! m3cor = soilfac*bv_soila + seasfac*bv_seas + & ! anthfac*bv_antha ! Cannot set_sulf here because it is a "radm2" species whose bc value ! is set via bdy_chem_value. Instead, xl(iref(p_sulf-1),:) is set to ! the value conmin in subroutine gasprofile_init_pnnl ! if( nch == p_sulf ) chem = conmin !as per rce's 0 recommendation END SUBROUTINE sorgam_set_aer_bc_pnnl !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE sorgam_vbs_set_aer_bc_pnnl( chem, z, nch, config_flags ) ! USE module_data_sorgam_vbs, ONLY : dginia, dginin, dginic, esn36, esc36, esa36, seasfac, no3fac, nh4fac, so4fac, soilfac, anthfac, orgfac implicit none INTEGER,INTENT(IN ) :: nch real,intent(in ) :: z REAL,INTENT(INOUT ) :: chem TYPE (grid_config_rec_type) , INTENT (in) :: config_flags REAL :: mult, & m3acc, m3cor, m3nuc, & bv_so4ai, bv_so4aj, & bv_nh4ai, bv_nh4aj, & bv_no3ai, bv_no3aj, & bv_eci, bv_ecj, & bv_p25i, bv_p25j, & bv_orgpai,bv_orgpaj, & bv_antha, bv_seas, bv_soila END SUBROUTINE sorgam_vbs_set_aer_bc_pnnl !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE module_input_smoke_data