MODULE aerosol_cloud_mod use fms_mod, only : error_mesg, FATAL, mpp_pe, & mpp_root_pe, open_namelist_file, & check_nml_error, close_file, & write_version_number, file_exist, & stdlog use constants_mod, ONLY : grav, cp_air, rdgas, rvgas, tfreeze use rad_utilities_mod, ONLY : aerosol_type use mpp_mod, only : mpp_clock_id, mpp_clock_begin, & mpp_clock_end, CLOCK_LOOP, & input_nml_file use polysvp_mod, ONLY : polysvp_l, polysvp_i, polysvp_init, & polysvp_end use aerosol_params_mod, only : aerosol_params_init, Nfact_du1, & Nfact_du2, Nfact_du3, Nfact_du4, & Nfact_du5, rbar_du1, rbar_du2, & rbar_du3, rbar_du4, rbar_du5 USE aer_ccn_act_mod, ONLY : aer_ccn_act_wpdf_m, & aer_ccn_act_init, aer_ccn_act_end USE ice_nucl_mod, ONLY : ice_nucl_wpdf, ice_nucl_wpdf_init, & ice_nucl_wpdf_end USE strat_cloud_utilities_mod, ONLY: strat_cloud_utilities_init, & diag_id_type, diag_pt_type, & strat_nml_type, atmos_state_type, & particles_type, strat_constants_type implicit none private !-------------------------------------------------------------------------- !---interfaces------------------------------------------------------------- public aerosol_cloud, aerosol_cloud_init, aerosol_cloud_end private aerosol_effects !-------------------------------------------------------------------------- !---version number--------------------------------------------------------- Character(len=128) :: Version = '$Id: aerosol_cloud.F90,v 20.0 2013/12/13 23:21:48 fms Exp $' Character(len=128) :: Tagname = '$Name: tikal $' !-------------------------------------------------------------------------- !---namelist--------------------------------------------------------------- integer :: rh_act_opt = 2 ! option as to which rh to use for ice ! nucleation. 1 => use gridbox mean; ! 2 => use rh in non-cloudy portion of ! box when the cloudy part of grid box ! is less than cf_thresh_nucl; use grid ! box mean otherwise. real :: sea_salt_scale_onl = 1. ! scaling factor to convert seasalt ! aerosol tracer to seasalt available ! for use as condensation nuclei logical :: reproduce_rk = .true. ! option to force use of same entries in ! aerosol activation tables as were used ! in legacy code; new code fixes an ! index offset error. real :: var_limit_ice = -999. ! lower limit to the vertical velocity ! variance of the pdf used for ice ! nucleation; default value implies ! that no limit is imposed. real :: cf_thresh_nucl = 0.98 ! threshold cloud fraction below which ! the rh of the non-cloudy portion of ! the gribox is used in determining ! ice nucleation (when rh_act_opt is >1) namelist / aerosol_cloud_nml / rh_act_opt, sea_salt_scale_onl, & reproduce_rk, var_limit_ice, & cf_thresh_nucl !------------------------------------------------------------------------- real, parameter :: d622 = rdgas / rvgas real, parameter :: d378 = 1. - d622 integer, parameter :: n_totmass = 4 integer, parameter :: n_imass = 12 LOGICAL :: do_var_limit_ice INTEGER :: wpdf_offs integer :: aero_init, aero_effects, aero_dust, aero_loop1, & aero_loop2, aero_loop3 logical :: module_is_initialized = .false. real :: missing_value = -1.e30 CONTAINS !######################################################################### subroutine aerosol_cloud_init (Constants) !------------------------------------------------------------------------ type(strat_constants_type), intent(in) :: Constants !-----local variables integer :: unit, io, ierr, logunit !----------------------------------------------------------------------- if (module_is_initialized) return !------------------------------------------------------------------------- ! process namelist. !------------------------------------------------------------------------- #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=aerosol_cloud_nml, iostat=io) ierr = check_nml_error(io,'aerosol_cloud_nml') #else if ( file_exist('input.nml')) then unit = open_namelist_file ( ) ierr=1; do while (ierr /= 0) read (unit, nml=aerosol_cloud_nml, iostat=io, end=10) ierr = check_nml_error(io,'aerosol_cloud_nml') enddo 10 call close_file (unit) endif #endif !------------------------------------------------------------------------- ! write version and namelist to standard log. !------------------------------------------------------------------------- call write_version_number(version, tagname) logunit = stdlog() if (mpp_pe() == mpp_root_pe()) & write (logunit, nml=aerosol_cloud_nml) if (var_limit_ice .EQ. -999. ) THEN do_var_limit_ice = .false. ELSE do_var_limit_ice = .true. END IF !------------------------------------------------------------------------- ! define offset into aerosol activation tables. !------------------------------------------------------------------------- IF (.NOT. (Constants%do_mg_microphys .OR. & Constants%do_mg_ncar_microphys) & .AND. reproduce_rk) THEN wpdf_offs =0 ELSE wpdf_offs = 1 END IF !------------------------------------------------------------------------- ! be sure needed modules have been initialized. !------------------------------------------------------------------------- call strat_cloud_utilities_init call polysvp_init call aer_ccn_act_init call ice_nucl_wpdf_init !------------------------------------------------------------------------- ! initialize clocks. !------------------------------------------------------------------------- aero_init = mpp_clock_id ('aerosol_cloud: initial', & grain = CLOCK_LOOP) aero_effects = mpp_clock_id ('aerosol_cloud: effects', & grain = CLOCK_LOOP) aero_dust = mpp_clock_id ('aerosol_cloud: dust ', & grain = CLOCK_LOOP) aero_loop1 = mpp_clock_id ('aerosol_cloud: loop1 ', & grain = CLOCK_LOOP) aero_loop2 = mpp_clock_id ('aerosol_cloud: loop2 ', & grain = CLOCK_LOOP) aero_loop3 = mpp_clock_id ('aerosol_cloud: loop3 ', & grain = CLOCK_LOOP) !------------------------------------------------------------------------- ! declare module initialized. !------------------------------------------------------------------------- module_is_initialized = .true. end subroutine aerosol_cloud_init !######################################################################### subroutine aerosol_cloud (idim, jdim, kdim, n_diag_4d, Nml, Constants, & Atmos_state, Particles, qa_upd, Aerosol, & diag_4d, diag_id, diag_pt, otun) !------------------------------------------------------------------------- integer, intent(in) :: idim, jdim, kdim, n_diag_4d type(strat_nml_type), intent(in) :: Nml type(strat_constants_type), intent(in) :: Constants type(atmos_state_type), intent(inout) :: Atmos_state type(particles_type), intent(inout) :: Particles REAL, DIMENSION(:,:,:), INTENT(IN ) :: qa_upd TYPE(aerosol_type), INTENT (in) :: Aerosol REAL, DIMENSION( :,:,:,0:), INTENT(INOUT ) :: diag_4d TYPE(diag_id_type), intent(in) :: diag_id TYPE(diag_pt_type), intent(in) :: diag_pt INTEGER, INTENT(IN ) :: otun !------------------------------------------------------------------------- !----local variables real, dimension(idim, jdim, kdim, n_totmass) :: totalmass1 real, dimension(idim, jdim, kdim, n_imass) :: imass1 real, dimension(idim, jdim, kdim) :: & up_strat, wp2, wp2i, wp2t, thickness, & ni_sulf, ni_dust, ni_bc, cf, u_i, u_l, & eslt, esit, qvsl, qvsi, qs_d, qvt INTEGER :: i, j, k !------------------------------------------------------------------------- ! initialize fields. !------------------------------------------------------------------------- call mpp_clock_begin (aero_init) totalmass1 = 0. imass1 = 0. call mpp_clock_end (aero_init) !--------------------------------------------------------------------- ! call aerosol_effects to determine aerosol fields that will impact the ! cloud / ice droplet numbers and the bergeron process, if these effects ! have been activated. !--------------------------------------------------------------------- call mpp_clock_begin (aero_effects) if (Nml%do_liq_num .or. Nml%do_dust_berg) then call aerosol_effects (idim, jdim, kdim, n_diag_4d, & Atmos_state%phalf, Atmos_state%airdens, & Atmos_state%T_in, Particles%concen_dust_sub,& totalmass1, imass1, Aerosol, Nml, diag_4d, & diag_id, diag_pt, otun ) endif call mpp_clock_end (aero_effects) !------------------------------------------------------------------------- ! define the number and mean radius of dust particles present. !------------------------------------------------------------------------- call mpp_clock_begin (aero_dust) do k = 1,kdim do j=1,jdim do i = 1,idim Particles%ndust(i,j,k) = Nfact_du1*imass1(i,j,k,8) + & Nfact_du2*imass1(i,j,k,9) + & Nfact_du3*imass1(i,j,k,10) + & Nfact_du4*imass1(i,j,k,11) + & Nfact_du5*imass1(i,j,k,12) Particles%rbar_dust(i,j,k) = ( & Nfact_du1*imass1(i,j,k,8)*rbar_du1 + & Nfact_du2*imass1(i,j,k,9)*rbar_du2 + & Nfact_du3*imass1(i,j,k,10)*rbar_du3 + & Nfact_du4*imass1(i,j,k,11)*rbar_du4 + & Nfact_du5*imass1(i,j,k,12)*rbar_du5 ) / & MAX (Particles%ndust(i,j,k),1.e-10) end do end do end do call mpp_clock_end (aero_dust) !------------------------------------------------------------------------- ! do the following calculations if droplet number is being predicted: !------------------------------------------------------------------------- if (Nml%do_liq_num) then !------------------------------------------------------------------------- ! compute the relevant upward velocity for droplet / ice activation. !------------------------------------------------------------------------- call mpp_clock_begin (aero_loop1) do k = 1,kdim do j=1,jdim do i = 1,idim up_strat(i,j,k) = -1.*(((Atmos_state%omega(i,j,k) + & grav*Atmos_state%Mc(i,j,k))/ & Atmos_state%airdens(i,j,k)/grav) + & Atmos_state%radturbten2(i,j,k)* & cp_air/grav) end do end do end do call mpp_clock_end (aero_loop1) !------------------------------------------------------------------------- ! define the layer thickness and variance of the vertical velocity. ! call aer_ccn_act_wpdf_m to determine number of activated droplets. ! this call need not be made if the activation vertical velocity is ! downward, the rotstayn-klein microphysics is being used, and pdf_clouds ! are not activated; in such a case, no particles are activated. !------------------------------------------------------------------------- call mpp_clock_begin (aero_loop2) do k = 1,kdim do j=1,jdim do i = 1,idim if ( (Nml%do_pdf_clouds) .or. & (Constants%do_mg_microphys) .or. & (Constants%do_mg_ncar_microphys) .or. & (up_strat(i,j,k) >= 0.0) ) then thickness(i,j,k) = Atmos_state%deltpg(i,j,k)/ & Atmos_state%airdens(i,j,k) wp2t(i,j,k) = 2.0/(3.0*0.548**2)* & (0.5*(Atmos_state%diff_t(i,j,k) + & Atmos_state%diff_t(i,j,min(k+1,KDIM)))/& thickness(i,j,k) )**2 wp2(i,j,k) = MAX (wp2t(i,j,k), Nml%var_limit**2) if(diag_id%subgrid_w_variance > 0) & diag_4d(i,j,k,diag_pt%subgrid_w_variance) = wp2(i,j,k)**0.5 call aer_ccn_act_wpdf_m & (Atmos_state%T_in(i,j,k), Atmos_state%pfull(i,j,k), & up_strat(i,j,k), wp2(i,j,k), wpdf_offs, & totalmass1(i,j,k,:), Particles%drop1(i,j,k)) else Particles%drop1(i,j,k) = 0. endif end do end do end do call mpp_clock_end (aero_loop2) !------------------------------------------------------------------------- ! if ice nuclei activation is being done via a vertical velocity pdf, ! execute the following code. !------------------------------------------------------------------------- call mpp_clock_begin (aero_loop3) if (Nml%do_ice_nucl_wpdf) THEN do k = 1,kdim do j = 1,jdim do i = 1,idim !------------------------------------------------------------------------ ! ice nuclei activation may only occur at temps below -5 C. !------------------------------------------------------------------------ if (Atmos_state%T_in(i,j,k) .LT. tfreeze - 5.) then !----------------------------------------------------------------------- ! place a lower limit on the velocity pdf variance, if desired. !----------------------------------------------------------------------- IF (do_var_limit_ice) THEN wp2i(i,j,k) = MAX (wp2t(i,j,k), var_limit_ice**2) ELSE wp2i(i,j,k) = wp2(i,j,k) END IF !------------------------------------------------------------------------ ! define the saturation specific humidities over liquid and ice. !------------------------------------------------------------------------ eslt(i,j,k) = polysvp_l(Atmos_state%T_in(i,j,k)) qs_d(i,j,k) = Atmos_state%pfull(i,j,k) - d378*eslt(i,j,k) qs_d(i,j,k) = max(qs_d(i,j,k),eslt(i,j,k)) qvsl(i,j,k) = 0.622 *eslt(i,j,k)/qs_d(i,j,k) esit(i,j,k) = polysvp_i(Atmos_state%T_in(i,j,k)) qs_d(i,j,k) = Atmos_state%pfull(i,j,k) - d378*esit(i,j,k) qs_d(i,j,k) = max(qs_d(i,j,k),esit(i,j,k)) qvsi(i,j,k) = 0.622 *esit(i,j,k)/qs_d(i,j,k) !------------------------------------------------------------------------ ! define the relative humidities wrt ice and liquid to be used in the ! calculation of ice nuclei activation. it may vary based on the nml ! variable rh_act_opt. !------------------------------------------------------------------------ IF (rh_act_opt .EQ. 1) THEN qvt(i,j,k) = Atmos_state%qv_In(i,j,k) cf(i,j,k) = 0. ELSE cf(i,j,k) = qa_upd(i,j,k)+ Atmos_state%ahuco(i,j,k) IF (cf(i,j,k) .LT. cf_thresh_nucl) THEN qvt(i,j,k) = (Atmos_state%qv_in(i,j,k) - & cf(i,j,k)*Atmos_State%qs(i,j,k))/ & (1. - cf(i,j,k)) ELSE qvt(i,j,k) = Atmos_state%qv_in(i,j,k) ENDIF END IF if (qvt(i,j,k) .LE. 0.0) then qvt(i,j,k) = MAX(Atmos_state%qv_in(i,j,k), Nml%qmin) endif u_i(i,j,k) = qvt(i,j,k)/qvsi(i,j,k) u_l(i,j,k) = qvt(i,j,k)/qvsl(i,j,k) !-------------------------------------------------------------------------- ! if debugging is active and the relative humidity exceeds 200%, output ! relevant variables. !-------------------------------------------------------------------------- if (Nml%debugo .and. & Atmos_state%T_in(i,j,k) .lt. 260. .and. & u_i(i,j,k) .gt. 200. ) then write(otun,*) " +++++++++++++++++++++++ " write(otun,*) " i,j,k ,u_i ", i,j,k,u_i(i,j,k) write(otun,*) " qs, qvsi , qvt, qv ", & Atmos_state%qs(i,j,k), qvsi(i,j,k), & qvt(i,j,k), Atmos_state%qv_in(i,j,k) write(otun,*) " cf, ahuco,qa_upd,qrat ", cf(i,j,k), & Atmos_state%ahuco(i,j,k), & qa_upd(i,j,k), Atmos_state%qrat(i,j,k) write(otun,*) " qv(i,k) - cf * qs(i,k), (1-cf) " , & Atmos_state%qv_in(i,j,k) - & cf(i,j,k)*Atmos_state%qs(i,j,k), & (1-cf(i,j,k)) write(otun,*) " +++++++++++++++++++++++ " endif !------------------------------------------------------------------------- ! call ice_nucl_wpdf to obtain number of activated ice crystals. !------------------------------------------------------------------------- call ice_nucl_wpdf ( & Atmos_state%T_in(i,j,k), u_i(i,j,k), u_l(i,j,k),& up_strat(i,j,k), wp2i(i,j,k), & Atmos_state%zfull(i,j,k), totalmass1(i,j,k,:), & imass1(i,j,k,:), n_totmass, n_imass, & Particles%crystal1(i,j,k), & Particles%drop1(i,j,k), Particles%hom(i,j,k), & Atmos_state%rh_crit(i,j,k), & Atmos_state%rh_crit_min(i,j,k), ni_sulf(i,j,k), & ni_dust(i,j,k), ni_bc(i,j,k)) else ni_sulf(i,j,k) = 0. ni_dust(i,j,k) = 0. ni_bc (i,j,k) = 0. cf(i,j,k) = missing_value u_i(i,j,k) = missing_value u_l(i,j,k) = missing_value endif !------------------------------------------------------------------------- ! define the critical relative humidity that was used for ice nuclei ! activation, averaged over all spectral intervals (rh_crit). also define ! the minimum value from any of the spectral regions (rh_crit_min). !------------------------------------------------------------------------- IF (Atmos_state%T_in(i,j,k) .LT. 250.) THEN Atmos_state%rh_crit(i,j,k) = & MAX(Atmos_state%rh_crit(i,j,k), 1.) Atmos_state%rh_crit_min(i,j,k) = & MAX(Atmos_state%rh_crit_min(i,j,k), 1.) !------------------------------------------------------------------------- ! if debugging is active, and the minimum is over 250%, output relevant ! information. !------------------------------------------------------------------------- IF ( Nml%debugo ) THEN IF (Atmos_state%rh_crit_min(i,j,k) .GT. 2.5 ) THEN write(otun,*) "MMMMMM RH", i,j,k write(otun,*) " rh_crit_min_1d ", & Atmos_state%rh_crit_min(i,j,k) END IF END IF ELSE !------------------------------------------------------------------------- ! if temp over 250 K, set values to 1.0. !------------------------------------------------------------------------- Atmos_state%rh_crit(i,j,k) = 1. Atmos_state%rh_crit_min(i,j,k) = 1. END IF end do end do end do !------------------------------------------------------------------------- ! define various desired diagnostics. !------------------------------------------------------------------------- if ( diag_id%imass7 > 0 ) & diag_4d(:,:,:,diag_pt%imass7) = imass1(:,:,:,7) if(diag_id%potential_crystals > 0) & diag_4d(:,:,:,diag_pt%potential_crystals) = & Particles%crystal1 if ( diag_id%rhcrit > 0 ) & diag_4d(:,:,:,diag_pt%rhcrit) = 100.*Atmos_state%rh_crit if ( diag_id%rhcrit_min > 0 ) & diag_4d(:,:,:,diag_pt%rhcrit_min) = 100.*Atmos_state%rh_crit_min if ( diag_id%ndust1 > 0 ) & diag_4d(:,:,:,diag_pt%ndust1) = Nfact_du1 * imass1(:,:,:,8) if ( diag_id%ndust2 > 0 ) & diag_4d(:,:,:,diag_pt%ndust2) = Nfact_du2 * imass1(:,:,:,9) if ( diag_id%ndust3 > 0 ) & diag_4d(:,:,:,diag_pt%ndust3) = Nfact_du3 * imass1(:,:,:,10) if ( diag_id%ndust4 > 0 ) & diag_4d(:,:,:,diag_pt%ndust4) = Nfact_du4 * imass1(:,:,:,11) if ( diag_id%ndust5 > 0 ) & diag_4d(:,:,:,diag_pt%ndust5) = Nfact_du5 * imass1(:,:,:,12) if ( diag_id%ni_dust > 0 ) & diag_4d(:,:,:,diag_pt%ni_dust) = ni_dust if ( diag_id%ni_bc > 0 ) diag_4d(:,:,:,diag_pt%ni_bc) = ni_bc if ( diag_id%ni_sulf > 0 ) & diag_4d(:,:,:,diag_pt%ni_sulf) = ni_sulf if (rh_act_opt .ne. 1) then if ( diag_id%cfin > 0 ) diag_4d(:,:,:,diag_pt%cfin) = cf end if if ( diag_id%rhiin > 0 ) diag_4d(:,:,:,diag_pt%rhiin) = u_i if ( diag_id%rhlin > 0 ) diag_4d(:,:,:,diag_pt%rhlin) = u_l END IF ! do_ice_nucl_wpdf if(diag_id%potential_droplets > 0) & diag_4d(:,:,:,diag_pt%potential_droplets) = & Particles%drop1 call mpp_clock_end (aero_loop3) end if ! (do_liq_num) !------------------------------------------------------------------------- END SUBROUTINE aerosol_cloud !##################################################################### subroutine aerosol_cloud_end module_is_initialized = .false. call polysvp_end call aer_ccn_act_end call ice_nucl_wpdf_end end subroutine aerosol_cloud_end !######################################################################### subroutine aerosol_effects (idim,jdim,kdim,n_diag_4d, phalf, airdens, T, & concen_dust_sub, totalmass1, imass1, Aerosol, & Nml, diag_4d, diag_id, diag_pt, otun ) INTEGER, INTENT (in) :: idim, jdim, & kdim, n_diag_4d REAL, DIMENSION(idim,jdim,kdim+1), INTENT(in ) :: phalf REAL, DIMENSION(idim,jdim,kdim), INTENT(in ) :: airdens, T REAL, DIMENSION(idim,jdim,kdim), INTENT(inout) :: concen_dust_sub REAL, DIMENSION(idim,jdim,kdim,n_totmass), INTENT(inout) :: totalmass1 REAL, DIMENSION(idim,jdim,kdim,n_imass), INTENT(inout) :: imass1 type(strat_nml_type), intent(in) :: Nml TYPE(aerosol_type), INTENT (in) :: Aerosol REAL, DIMENSION(idim,jdim,kdim,0:n_diag_4d), & INTENT(INOUT) :: diag_4d TYPE(diag_id_type), intent(in) :: diag_id TYPE(diag_pt_type), intent(in) :: diag_pt INTEGER, INTENT (in) :: otun REAL, DIMENSION(idim,jdim,kdim) :: pthickness, concen_all_sub, & concen_ss_sub, concen_ss_sup INTEGER :: i, j, k, na, s LOGICAL :: used !------------------------------------------------------------------------- ! initialize output field (sub-micron dust content). !------------------------------------------------------------------------- do k=1,kdim do j=1,jdim do i=1,idim concen_dust_sub(i,j,k) = 0. end do end do end do !------------------------------------------------------------------------- ! initialize local accumulation arrays. !------------------------------------------------------------------------- if (Nml%use_online_aerosol) then do k=1,kdim do j=1,jdim do i=1,idim concen_ss_sub(i,j,k) = 0. concen_ss_sup(i,j,k) = 0. concen_all_sub(i,j,k) = 0. end do end do end do endif !------------------------------------------------------------------------ ! define layer thickness. !------------------------------------------------------------------------ do k=1,kdim do j=1,jdim do i=1,idim if (phalf(i,j,k) < 1.0) then pthickness(i,j,k) = (phalf(i,j,k+1) - phalf(i,j,k))/& grav/airdens(i,j,k) else pthickness(i,j,k) = log(phalf(i,j,k+1)/ & phalf(i,j,k))*8.314*T(i,j,k)/(9.8*0.02888) end if end do end do end do !------------------------------------------------------------------------- ! define the aerosol content in various categories that may be used for ! droplet or ice crystal nucleation. save various diagnostics as they ! are calculated. !------------------------------------------------------------------------- if (Nml%do_liq_num) then if (Nml%use_online_aerosol) then do na = 1,size(Aerosol%aerosol,4) if (trim(Aerosol%aerosol_names(na)) == 'so4' .or. & trim(Aerosol%aerosol_names(na)) == 'so4_anthro' .or.& trim(Aerosol%aerosol_names(na)) == 'so4_natural') then do k=1,kdim do j=1,jdim do i=1,idim totalmass1(i,j,k,1) = totalmass1(i,j,k,1) + & Aerosol%aerosol(i,j,k,na) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'omphilic' .or.& trim(Aerosol%aerosol_names(na)) == 'omphobic') then do k=1,kdim do j=1,jdim do i=1,idim totalmass1(i,j,k,4) = totalmass1(i,j,k,4) + & Aerosol%aerosol(i,j,k,na) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'seasalt1' ) then do k=1,kdim do j=1,jdim do i=1,idim concen_ss_sub(i,j,k) = concen_ss_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,1) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'seasalt2') then do k=1,kdim do j=1,jdim do i=1,idim concen_ss_sub(i,j,k) = concen_ss_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,2) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'seasalt3') then do k=1,kdim do j=1,jdim do i=1,idim concen_ss_sup(i,j,k) = concen_ss_sup(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,3) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'seasalt4' ) then do k=1,kdim do j=1,jdim do i=1,idim concen_ss_sup(i,j,k) = concen_ss_sup(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,4) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'seasalt5' ) then do k=1,kdim do j=1,jdim do i=1,idim concen_ss_sup(i,j,k) = concen_ss_sup(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,5) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'bcphilic') then do k=1,kdim do j=1,jdim do i=1,idim concen_all_sub(i,j,k) = concen_all_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,6) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'bcphobic' ) then do k=1,kdim do j=1,jdim do i=1,idim concen_all_sub(i,j,k) = concen_all_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'dust1') then do k=1,kdim do j=1,jdim do i=1,idim concen_all_sub(i,j,k) = concen_all_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,8) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) ! according to Paul dust 1 and dust2 are sub- 2.5 micron imass1(i,j,k,7) = imass1(i,j,k,7) + & Aerosol%aerosol(i,j,k,na) if (Nml%do_dust_berg) then concen_dust_sub(i,j,k) = concen_dust_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) endif end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'dust2' ) then do k=1,kdim do j=1,jdim do i=1,idim concen_all_sub(i,j,k) = concen_all_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,9) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) ! according to Paul dust 1 and dust2 are sub- 2.5 micron imass1(i,j,k,7) = imass1(i,j,k,7) + & Aerosol%aerosol(i,j,k,na) if (Nml%do_dust_berg) then concen_dust_sub(i,j,k) = concen_dust_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) endif end do end do end do else if(trim(Aerosol%aerosol_names(na)) == 'dust3' ) then do k=1,kdim do j=1,jdim do i=1,idim concen_all_sub(i,j,k) = concen_all_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) imass1(i,j,k,10) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) if (Nml%do_dust_berg) then concen_dust_sub(i,j,k) = concen_dust_sub(i,j,k) + & Aerosol%aerosol(i,j,k,na) endif end do end do end do else if (trim(Aerosol%aerosol_names(na)) == 'dust4') then do k=1,kdim do j=1,jdim do i=1,idim imass1(i,j,k,11) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) end do end do end do else if (trim(Aerosol%aerosol_names(na)) == 'dust5') then do k=1,kdim do j=1,jdim do i=1,idim imass1(i,j,k,12) = Aerosol%aerosol(i,j,k,na)/ & pthickness(i,j,k) end do end do end do endif end do do k=1,kdim do j=1,jdim do i=1,idim totalmass1(i,j,k,3) = concen_ss_sub(i,j,k) totalmass1(i,j,k,2) = concen_all_sub(i,j,k) + & totalmass1(i,j,k,4) + & concen_ss_sub(i,j,k) if (Nml%use_sub_seasalt) then else totalmass1(i,j,k,3) = concen_ss_sub(i,j,k) + & concen_ss_sup(i,j,k) endif if (sea_salt_scale_onl .NE. 1. ) then totalmass1(i,j,k,3) = sea_salt_scale_onl* & totalmass1(i,j,k,3) end if if (diag_id%sulfate > 0) then diag_4d(i,j,k, diag_pt%sulfate) = & 0.7273*totalmass1(i,j,k,1)/ & pthickness(i,j,k)*1.0e9 endif if (diag_id%seasalt_sub > 0) then diag_4d(i,j,k,diag_pt%seasalt_sub) = & concen_ss_sub(i,j,k)/ & pthickness(i,j,k)*1.0e9 endif if (diag_id%seasalt_sup > 0) then diag_4d(i,j,k, diag_pt%seasalt_sup) = & concen_ss_sup(i,j,k)/ & pthickness(i,j,k)*1.0e9 endif end do end do end do else ! (use_online_aerosol) if (Nml%do_dust_berg) then ! YMice submicron dust (NO. 14 to NO. 18) do s = 14,18 do k=1,kdim do j=1,jdim do i=1,idim concen_dust_sub(i,j,k) = concen_dust_sub(i,j,k)+ & Aerosol%aerosol(i,j,k,s) end do end do end do end do endif if (diag_id%sulfate > 0) then ! anthro. and natural sulfate concentration (ug so4/m3) do k=1,kdim do j=1,jdim do i=1,idim diag_4d(i,j,k, diag_pt%sulfate) = 0.7273* & (Aerosol%aerosol(i,j,k,1) + & Aerosol%aerosol(i,j,k,2))/& pthickness(i,j,k)*1.0e9 end do end do end do endif !offline ! NO. 1 natural Sulfate; NO. 2 anthro. sulfate; NO. 3 Sea Salt; NO. 4 Or ganics do k=1,kdim do j=1,jdim do i=1,idim totalmass1(i,j,k,1) = Aerosol%aerosol(i,j,k,2) totalmass1(i,j,k,2) = Aerosol%aerosol(i,j,k,1) totalmass1(i,j,k,3) = Nml%sea_salt_scale* & Aerosol%aerosol(i,j,k,5) totalmass1(i,j,k,4) = Nml%om_to_oc* & Aerosol%aerosol(i,j,k,3) end do end do end do endif ! (use_online_aerosol) !cms2008-10-08 do na = 1, 3 do na = 1, 4 do k=1,kdim do j=1,jdim do i=1,idim totalmass1(i,j,k,na) = totalmass1(i,j,k,na)/ & pthickness(i,j,k)*1.0e9*1.0e-12 end do end do end do end do do k=1,kdim do j=1,jdim do i=1,idim ! submicron dust concentration (ug/m3) (NO. 2 to NO. 4) imass1(i,j,k,7) = imass1(i,j,k,7)/pthickness(i,j,k)*1.0e9 end do end do end do if (Nml%do_dust_berg) then ! submicron dust concentration (ug/m3) (NO. 2 to NO. 4) do k=1,kdim do j=1,jdim do i=1,idim concen_dust_sub(i,j,k) = concen_dust_sub(i,j,k)/ & pthickness(i,j,k)*1.0e9 end do end do end do endif if (diag_id%om > 0) then do k=1,kdim do j=1,jdim do i=1,idim !RSH the 1.0e12 factor here is to counter the 1.0e-12 factor applied ! above to totalmass1. The 1.0e9 factor (kg -> ug) is already in totalmass1. diag_4d(i,j,k,diag_pt%om) = totalmass1(i,j,k,2)*1.0e12 end do end do end do endif endif ! (do_liq_num) !---------------------------------------------------------------------- end subroutine aerosol_effects !######################################################################### END MODULE aerosol_cloud_mod