! ! Contact: thomas.diehl@nasa.gov ! ! Fortran subroutine stub for reading the volcanic SO2 emission data into a ! CTM or climate model. ! This routine assumes that the variables ! -> nevents, jdn, elevation, cloud_column_height, so2 <- ! from the volcanic emission files have been read into the variables ! -> nve, vjdn, velev, vcch, vso2 <- ! elsewhere. It also assumes that the grid indeces for each volcano have been ! computed from the variables lon and lat in the emission files elsewhere and ! stored in the variables vi and vj. ! It is also assumed that vso2 has been converted from kt SO2/event/day to ! kg SO2/event/day. ! ! imx: number of grid boxes in east-west direction ! jmx: number of grid boxes in north-south direction ! lmx: number of vertical grid boxes ! delz(imx,jmx,lmx): height of grid boxes [m] ! surface_altitude(imx,jmx): surface altitude [m] ! eso2_volc(imx,jmx,lmx): amount of so2 injected into a grid box [kg SO2/box/s] ! Module variables: ! LOGICAL :: levolc ! (volcanic emissions are computed if set to true) ! INTEGER :: idx1, next_idx1 ! LOGICAL :: first = .true. ! INTEGER :: nve ! allocate these variables with dimension nve: ! INTEGER, ALLOCATABLE :: vi(:), vj(:), vjdn(:), vcch(:), velev(:) ! REAL, ALLOCATABLE :: vso2(:), vlon(:), vlat(:) ! dt: dt(1)%jday contains the Julian Day Number at the beginning of the current ! model time step. ! jday_m1: Julian Day Number at beginning of previous time step SUBROUTINE src_vso2_stub ( imx, jmx, lmx, & delz, surface_altitude, & eso2_volc ) USE mo_control, ONLY: levolc USE mo_volc, ONLY: nve, vjdn, vcch, velev, vi, vj, vso2, & idx1, first, next_idx1 USE mo_time_control, ONLY: dt, jday_m1 USE mo_exception, ONLY: finish IMPLICIT NONE INTEGER, INTENT(IN) :: imx, jmx, lmx REAL, INTENT(IN) :: delz(imx,jmx,lmx) REAL, INTENT(IN) :: surface_altitude(imx,jmx) REAL, INTENT(OUT) :: eso2_volc(imx,jmx,lmx) INTEGER :: k, l, lv1, lv2, nvolc INTEGER :: idx1_tmp(1) REAL :: hght, slab, slab1 REAL :: zh(0:lmx), dz(lmx) LOGICAL :: found CHARACTER(LEN=100) :: msg ! executable statements ! **************************************************************************** ! * Calculate volcanic emission of SO2. * ! **************************************************************************** eso2_volc(:,:,:) = 0.0 if_volc_emis: IF (levolc) THEN IF (first) THEN nvolc = COUNT(vjdn == dt(1)%jday) IF (nvolc == 0) THEN msg = 'Current date not found in volcanic database.' CALL finish('source_su', msg) END IF idx1_tmp = MAXLOC(vjdn, MASK = vjdn == dt(1)%jday) idx1 = idx1_tmp(1) END IF IF (dt(1)%jday /= jday_m1) THEN idx1 = next_idx1 jday_m1 = dt(1)%jday END IF first = .FALSE. found = .FALSE. nevents_loop: DO k = idx1, nevents ! all volcanic events of one day are stored sequentially; therefore, ! once we encounter the first element which belongs to the next day ! the loop can be terminated IF (dt(1)%jday /= vjdn(k)) THEN next_idx1 = k EXIT END IF found = .TRUE. IF (velev(k) == vcch(k)) THEN ! the volcano inventory is constructed with vcch=velev for ! non-eruptive degassing (no plume is assumed); ! therefore, the emission is placed only in the level which ! contains the crater elevation hght = REAL(velev(k)) zh(0) = 0.0 + surface_altitude(vi(k),vj(k)) ! convert to m a.s.l. IF (hght < zh(0)) THEN eso2_volc(vi(k),vj(k),1) = & eso2_volc(vi(k),vj(k),1) + vso2(k)/(24.0*3600.0) ELSE DO l = 1,lmx dz(l) = delz(vi(k),vj(k),lmx-l+1) ! delz: TOA->sfc zh(l) = zh(l-1) + dz(l) IF (zh(l-1) <= hght .AND. zh(l) > hght) THEN eso2_volc(vi(k),vj(k),l) = eso2_volc(vi(k),vj(k),l) + & vso2(k)/(24.0*3600.0) EXIT END IF END DO END IF ELSE ! -------------------------------------------------- ! Define a slab at the top 1/3 of the volcano plume. ! -------------------------------------------------- hght = REAL(vcch(k)) slab1 = hght - (hght - REAL(velev(k)))/3.0 ! slab bottom height slab = hght - slab1 ! slab thickness zh(0) = 0.0 + surface_altitude(vi(k),vj(k)) ! convert to m a.s.l. IF (hght <= zh(0)) THEN lv1 = 1 lv2 = 1 dz(1) = slab ELSE level_loop: DO l = 1,lmx dz(l) = delz(vi(k),vj(k),lmx-l+1) ! delz: TOA->sfc zh(l) = zh(l-1) + dz(l) ! top edge height IF (l == lmx .AND. hght > zh(l)) THEN ! max model erup.height hght = zh(l) slab1 = slab1 - (hght - zh(l)) END IF ! -------------------------------------------------- ! If Zh(l) <= bottom of the slab, go to next level. ! -------------------------------------------------- IF (zh(l) <= slab1) CYCLE ! -------------------------------------------------- ! If the slab is only in current level: ! -------------------------------------------------- IF (zh(l-1) <= slab1 .AND. zh(l) >= hght) THEN lv1 = l lv2 = l dz(l) = slab EXIT ! -------------------------------------------------- ! The slab extends more than one level. Find the ! lowest (lv1) and the highest (lv2) levels: ! -------------------------------------------------- ELSE IF (zh(l-1) < slab1 .AND. zh(l) > slab1) THEN lv1 = l dz(l) = zh(l) - slab1 ELSE IF (zh(l-1) <= hght .AND. zh(l) > hght) THEN lv2 = l dz(l) = hght - zh(l-1) EXIT END IF END DO level_loop END IF ! -------------------------------------------------- ! Calculate SO2 emission in the levels between ! lv1 and lv2. Convert vso2 from kgSO2/box/day ! to kgSO2/box/sec. eso2_volc is distributed evenly ! with altitude among the slab. ! Summation is necessary since more than one volcano ! could be located within the same grid box. ! -------------------------------------------------- DO l = lv1,lv2 eso2_volc(vi(k),vj(k),l) = eso2_volc(vi(k),vj(k),l) & + vso2(k)/(24.0*3600.0) * dz(l)/slab END DO END IF END DO nevents_loop IF (.NOT. found) THEN msg = 'Current date not found in volcanic database.' CALL finish('source_su', msg) END IF END IF if_volc_emis END SUBROUTINE src_vso2_stub