#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: OPTICS ! ! !DESCRIPTION: Optics module to calculate optical depth from m7 output, ! based on the AOP_Package op Michael Kahnert. ! ! The optics code should be used in the following way (see ! examples in photolysis.F90, ecearth_optics.F90, and ! user_output_aerocom.F90): ! ! 1) prepare the optics calculation by calling ! 'OPTICS_INIT' --> lookuptables etc. ! this routine reads the wavelengths, lookupable and calculates ! refr. indices at these wavelengths. ! ! 2) allocate AOP arrays aop_out_ext/a/g with ! (n_gridcells, n_wavelengths, n_split), with: ! 'n_split' = 1 for (split == .false.) or ! 'n_split' = 11 for (split == .true. ), ie. ! partial contributions by ! Total, SO4, BC, OC, SS, DU, NO3, Water, Fine, Fine Dust, Fine SS ! additional fields are also provided for split==.true. in ! aop_out_add, which has to have size ! (n_gridcells, n_wavelengths, n_add) with nadd = 2 for ! surface PM10 dry extinction and surface PM10 dry absorption ! 3) Compute AOP for current conditions by calling 'OPTICS_AOP_GET' ! 4) done: 'OPTICS_DONE' ! ! IMPORTANT: *) Skipped the ZOOM options! (temporary) ! *) OC is actually POM. Remember converting OC to POM ! when sending it to optics_calculate_aop !\\ !\\ ! !INTERFACE: ! MODULE OPTICS ! ! !USES: ! use GO, only : gol, goErr, goPr Use mo_aero_m7, only : nmod Use global_data, only : rcfile use GO, only : TrcFile, Init, Done, ReadRc Use Dims, only : nregions IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: optics_init, optics_done ! init/done methods for a wavelendep object public :: optics_aop_get ! compute aerosol optical properties ! ! !PRIVATE DATA MEMBERS: ! type(TrcFile) :: rcF character(len=256) :: lookuptable, refractive_indices ! ! !PUBLIC TYPES: ! ! wavelength type (to be used by methods using the optics) type, public :: wavelendep real :: wl ! user requested wavelength unit = um (e.g. 0.550) real, dimension(7) :: n ! SO4, BC, OC, SOA, SS, DU, WATER real, dimension(7) :: k ! SO4, BC, OC, SOA, SS, DU, WATER logical :: split, insitu end type wavelendep ! ! !PRIVATE TYPES: ! ! AOP input type and field type aopi real, dimension(nmod) :: SO4, NO3, BC, OC, SOA, SS, DU, H2O, numdens, rg, rgd end type aopi type(aopi), dimension(:), allocatable :: aop_in ! characteristics of the lookup-tables integer, parameter :: n_rii = 15 integer, parameter :: n_rir = 40 integer, parameter :: n_x = 100 real,dimension(n_rii) :: lkval ! -log img part refr. index real,dimension(n_rii) :: kval ! img part refr. index, 10^(-lkval) real,dimension(n_rir) :: n1r ! real part refr. index Real, Dimension(n_x) :: xs Real, Dimension(n_x,n_rir,n_rii),Target :: cext_159, a_159, g_159 real, dimension(n_x,n_rir,n_rii),Target :: cext_200, a_200, g_200 character(len=*), parameter :: mname = 'Optics' ! lookup table information for interpolation to wavelength properties Integer, Parameter :: opacdim = 61 ! Wavelength dimension of OPAC data Integer, Parameter :: echamhamdim = 49 ! Wavelength dimension of ECHAM-HAM data Integer, Parameter :: segelsteindim = 1261 ! Wavelength dimension of Segelstein data Real, Dimension(:,:), allocatable :: opac, echamham, segelstein ! ! !REVISION HISTORY: ! Implemented by Maarten Krol 03/2009 ! Extended by Joost Aan de Brugh ! ! 6 Feb 2011 - Achim Strunk - complete revision ! - worked on DECOUPLING optics routines ! from (optics)_output, in order to ! (re)establish a flexible way of using it ! - remaining parts have been moved to ! (the new routine) user_output_optics ! ! 24 Jun 2011 - Achim Strunk - added NO3 explicitly in order to allow ! for a slightly better split of the ! optical properties of (SO4+NO3) ! 20 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OPTICS_INIT ! ! !DESCRIPTION: Initialize a wavelendep object: Input/lookuptables are ! opened and read in, and wavelength dependent parameters are ! calculated. !\\ !\\ ! !INTERFACE: ! SUBROUTINE OPTICS_INIT( nwav, wdep, status ) ! ! !USES: ! USE MDF, ONLY : MDF_OPEN, MDF_NETCDF, MDF_READ, MDF_Get_Var, MDF_Close, MDF_Inq_VarID ! ! !INPUT PARAMETERS: ! integer, intent(in) :: nwav ! ! !OUTPUT PARAMETERS: ! type(wavelendep), dimension(nwav), intent(inout) :: wdep integer, intent(out) :: status ! ! !REMARKS: read by each core (FIXME?) ! !EOP !------------------------------------------------------------------------ !BOC integer :: fid, varid, i character(len=*), parameter :: rname = mname//'/Optics_Init' ! --- begin -------------------------------- allocate( opac (7, opacdim ) ) allocate( echamham (5, echamhamdim ) ) allocate( segelstein(3, segelsteindim) ) ! get filenames from rcfile call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'optics.lookuptable', lookuptable, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'optics.refractive_indices', refractive_indices, status ) IF_NOTOK_RETURN(status=1) ! read lookup table CALL MDF_Open( TRIM(lookuptable), MDF_NETCDF, MDF_READ, fid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( fid, 'mr', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, n1r, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( fid, 'mi', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, kval, status ) IF_NOTOK_RETURN(status=1) lkval = -log10(kval) CALL MDF_Inq_VarID( fid, 'x', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, xs, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( fid, 'ext_159', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, cext_159, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( fid, 'ext_200', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, cext_200, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( fid, 'a_159', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, a_159, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( fid, 'a_200', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, a_200, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( fid, 'g_159', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, g_159, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Inq_VarID( fid, 'g_200', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, g_200, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Close( fid, status ) IF_NOTOK_RETURN(status=1) ! Read refractive indices CALL MDF_Open( TRIM(refractive_indices), MDF_NETCDF, MDF_READ, fid, status ) IF_NOTOK_RETURN(status=1) ! Get Opac data CALL MDF_Inq_VarID( fid, 'Opac', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, opac, status ) IF_NOTOK_RETURN(status=1) ! Get ECHAM-HAM data CALL MDF_Inq_VarID( fid, 'ECHAM-HAM', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, echamham, status ) IF_NOTOK_RETURN(status=1) ! Get Segelstein data CALL MDF_Inq_VarID( fid, 'Segelstein', varid, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( fid, varid, segelstein, status ) IF_NOTOK_RETURN(status=1) CALL MDF_Close( fid, status ) IF_NOTOK_RETURN(status=1) call optics_wavelen_init( wdep, status ) IF_ERROR_RETURN( status=1 ) write(gol,*) 'Optical parameters:'; call goPr do i = 1, nwav write(gol,*) 'Wavelength :', wdep(i)%wl ; call goPr write(gol,*) ' SO4 (real/img) :', wdep(i)%n(1), wdep(i)%k(1) ; call goPr write(gol,*) ' BC (real/img) :', wdep(i)%n(2), wdep(i)%k(2) ; call goPr write(gol,*) ' OC (real/img) :', wdep(i)%n(3), wdep(i)%k(3) ; call goPr write(gol,*) ' SOA (real/img) :', wdep(i)%n(4), wdep(i)%k(4) ; call goPr write(gol,*) ' SS (real/img) :', wdep(i)%n(5), wdep(i)%k(5) ; call goPr write(gol,*) ' DU (real/img) :', wdep(i)%n(6), wdep(i)%k(6) ; call goPr write(gol,*) ' H2O (real/img) :', wdep(i)%n(7), wdep(i)%k(7) ; call goPr enddo ! Done ! ------------- call Done( rcF, status ) IF_NOTOK_RETURN(status=1) deallocate( opac, echamham, segelstein ) status = 0 END SUBROUTINE OPTICS_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OPTICS_DONE ! ! !DESCRIPTION: dummy for now !\\ !\\ ! !INTERFACE: ! SUBROUTINE OPTICS_DONE( status ) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 6 Feb 2011 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC ! ok status = 0 END SUBROUTINE OPTICS_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OPTICS_WAVELEN_INIT ! ! !DESCRIPTION: Initialise parameters which are depending on given ! wavelengths. !\\ !\\ ! !INTERFACE: ! SUBROUTINE OPTICS_WAVELEN_INIT( wdep, status ) ! ! !INPUT/OUTPUT PARAMETERS: ! ! wavelength properties (wavelength itself and real/img part of refractive index) type(wavelendep), intent(inout), dimension(:) :: wdep ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 6 Feb 2011 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/optics_wavelen_init' integer :: nwl, i, idx real :: wl, h real :: nscale, kscale ! Real part n and imaginary part k of refractive index ! for each wavelength: !>>> TvN ! The refractive index for 'sulfate' is taken from ! the OPAC database (Hess et al., 1998; ! Koepke et al., MPI report no. 243, 1997). ! The value used is that of 'sulfate solution', ! i.e. particles consisting of 75% of sulfuric acid (H2SO4), ! based on Fenn et al. (chapter 18 in Handbook of Geophysics ! and Space Environment, 1985). ! This is in line with Kinne et al. (JGR, 2003), ! who write that refractive indices for sulfate ! are usually based on 75% sulfuric acid solution. ! Actually, the OPAC refractive index agrees very well ! with the expression given by Boucher and Anderson (JGR, 1995) ! for H2SO4 at visible wavelengths. ! Thus, the OPAC data can be considered to apply to pure H2SO4, ! which is consistent of our application of mixing rules. ! For BC, the OPAC refractive index is scaled to one of ! the values proposed by Bond and Bergstrom (2006), ! valid at 550 nm (see their Table 5). ! Selected values for high, medium and low absorption are ! 1.95 + 0.79 i, 1.85 + 0.71 i, and 1.75 + 0.63 i, respectively. ! We select the low-absorption value, ! because it produces results in best agreement ! with AAOD from MODIS Collection 6 Deep Blue. ! Note that the medium-absorption value 1.85 + 0.71 i ! was used in ECHAM simulations by Stier et al. (2007), ! and found to give reasonable results. ! Lowenthal et al. (Atmos. Environ., 2000) ! give a value of 1.90 + 0.6 i. ! The reference value used in the scaling ! is the OPAC value at 550 nm: 1.75 + 0.44 i. ! It would be better to include the scaling ! in the input file. !nscale = 1.0 ! for using OPAC !kscale = 1.0 ! for using OPAC !nscale = 1.95/1.75 !kscale = 0.79/0.44 nscale = 1.85/1.75 kscale = 0.71/0.44 !nscale = 1.75/1.75 !kscale = 0.63/0.44 !<<< TvN nwl = size(wdep) do i = 1, nwl wl = wdep(i)%wl ! Interpolate Opac data findOpac: Do idx = 1,opacdim If(opac(1,idx) .gt. wl) exit findOpac End Do findOpac If(idx .gt. opacdim) then idx = opacdim h = 1.0 Else If(idx .eq. 1) then idx = 2 h = 0.0 Else h = (wl-opac(1,idx-1))/(opac(1,idx)-opac(1,idx-1)) End If ! SO4 wdep(i)%n(1) = opac(2,idx-1)+h*(opac(2,idx)-opac(2,idx-1)) wdep(i)%k(1) = opac(3,idx-1)+h*(opac(3,idx)-opac(3,idx-1)) ! BC !>>> TvN !wdep(i)%n(2) = opac(4,idx-1)+h*(opac(4,idx)-opac(4,idx-1)) !wdep(i)%k(2) = opac(5,idx-1)+h*(opac(5,idx)-opac(5,idx-1)) wdep(i)%n(2) = nscale * ( opac(4,idx-1)+h*(opac(4,idx)-opac(4,idx-1)) ) wdep(i)%k(2) = kscale * ( opac(5,idx-1)+h*(opac(5,idx)-opac(5,idx-1)) ) !<<< TvN ! SS wdep(i)%n(5) = opac(6,idx-1)+h*(opac(6,idx)-opac(6,idx-1)) wdep(i)%k(5) = opac(7,idx-1)+h*(opac(7,idx)-opac(7,idx-1)) ! Interpolate ECHAM-HAM data findEchamham: Do idx = 1,echamhamdim If(echamham(1,idx) .gt. wl) exit findEchamham End Do findEchamham If(idx .gt. echamhamdim) then idx = echamhamdim h = 1.0 Else If(idx .eq. 1) then idx = 2 h = 0.0 Else h = (wl-echamham(1,idx-1))/(echamham(1,idx)-echamham(1,idx-1)) End If ! OC !>>> TvN ! The 'ECHAM-HAM' data currently used for POM ! are based on the data from Fenn et al. (1985) ! for the 'water soluble' component, ! but at a reduced number of wavelengths up to 15 um. ! It would be better to use the original table ! in the input file. !<<< TvN wdep(i)%n(3) = echamham(2,idx-1)+h*(echamham(2,idx)-echamham(2,idx-1)) wdep(i)%k(3) = echamham(3,idx-1)+h*(echamham(3,idx)-echamham(3,idx-1)) ! SOA ! >>TB ! For now (March 2017) we use OC indices for SOA ! < 1 ) then if ( abs(n-n1r(i_n-1)) < abs(n-n1r(i_n)) ) i_n = i_n - 1 else if ( i_n < n_rir ) then if ( abs(n-n1r(i_n+1)) < abs(n-n1r(i_n)) ) i_n = i_n + 1 end if do i_k = 1, n_rii if (abs(lk - lkval(i_k)) < 1e-4) exit enddo ! ! PLS - test : ik can be determined without the preceding loop "do i_k = 1, n_rii" ! if (i_k.ne.i_knew) then ! status = 1 ! write (gol,'(" PLSPLS ik NE iknew = ",2(i2,2x))')i_k,i_knew ! endif ! following the "get_k" loop above, the only way to get into this is to ! have a NaN for k in the first place ? if (i_k > n_rii) then status = 1 write(gol,*)'FATAL ERROR: i_k value outside LUT' ! write (gol,'(" lk = ",E16.4)')lk ! write (gol,'(" k1 = ",E16.4)')k1 ! write (gol,'(" k = ",E16.4)')k ! write (gol,'(" n = ",E16.4)')n ! do i_n = 1, 15 ! write (gol,'(" lkval(",i2,") = ",E16.4)')i_n,lkval(i_n) ! call goPr ! write (gol,'(" kval(",i2,") = ",E16.4)')i_n,kval(i_n) ! call goPr ! enddo return endif ! Added check (15-7-2010 - P. Le Sager) if (i_n > n_rir) then status = 1 write(gol,*)'FATAL ERROR: i_n value out of range' return endif !>>> TvN ! Since xs(1) equals zero a problem may occur when xg is negative, ! because modrad.gt.xg would become true for i=1 in that case, ! while modrad1, Cext1, a1 and g1 are not yet set. ! It is not clear if negative xg values can ever occur, ! but if they do that should be prevented when calculating rg. ! ! However, the problem may perhaps already occur ! when xg equals zero, because of the finite machine precision. ! ! In any case, it is desired to initialize modrad1, Cext1, a1 and g1. ! Cext1, a1 and g1 should be set to their table entries for i=1, ! which are all zero: Cext1 = cext_table(1,i_n,i_k) a1 = a_table(1,i_n,i_k) g1 = g_table(1,i_n,i_k) ! modrad1 can be initialized to any value different from xs(1), ! to prevent division by zero. modrad1=xs(1)-9.99e-4 ! This combination will force slope to zero ! and Cext, a and g to the table entries for i=1 (zero), ! in case modrad.gt.xg is true for i=1. !<<< TvN get_values: do i = 1, n_x modrad = xs(i) a = a_table(i,i_n, i_k) g = g_table(i,i_n, i_k) cext = cext_table(i,i_n, i_k) ! With small concentrations, like in the stratosphere, M7 does not trust its radii. ! See m7_averageproperties -> zinsvol. The M7-radius goes to zero ! Modrad may never be larger than xg the first time. Modrad1 is not set, ! will be zero (or something worse) and dr will be zero (or something worse). ! slope is demolished, makeing all values NaN. Therefore, it is modrad .gt. xg ! instead of modrad .ge. xg ! ! PLS-TRANSLATION - It boils down to: "on the first iteration of this loop, if modrad=xg and ! you go into the if-statement below, you are in trouble, because modrad1 can be ! anything. To prevent that, replace GE with GT to avoid going into the if-statement at the ! first iteration." ! if(modrad.gt.xg)then dr=modrad-modrad1 dr1=xg-modrad1 slope=(Cext-Cext1)/dr Cext=Cext1+slope*dr1 slope=(a-a1)/dr a=a1+slope*dr1 slope=(g-g1)/dr g=g1+slope*dr1 exit get_values endif modrad1=modrad Cext1=Cext a1=a g1=g enddo get_values END SUBROUTINE OPTICS_GET !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: GET_REFR_IDX ! ! !DESCRIPTION: Compute refractive index of internally mixed aerosols by use ! of effective medium theory for the size-dependent aerosol ! mixtures assumed in M7. !\\ !\\ ! !INTERFACE: ! PURE SUBROUTINE GET_REFR_IDX(wdep, SO4, BC, OC, SOA, SS, DU, water, mode, m_eff, status, gol) ! ! !USES: ! use binas, only: rol use chem_param, only: ss_density, dust_density, carbon_density use chem_param, only: pom_density, so4_density, soa_density ! ! !INPUT PARAMETERS: ! type(wavelendep), intent(in) :: wdep ! wavelength properties (wavelength, re/img part of refractive index) real, intent(in) :: SO4, BC, OC,SOA ! mass mixing ratios or concentrations of sulphate, black carbon, organic carbon real, intent(in) :: SS, DU, water ! sea salt, dust, and water integer, intent(in) :: mode ! mode number (M7) ! ! !OUTPUT PARAMETERS: ! complex, intent(out) :: m_eff ! effective refractive index of mixture integer, intent(out) :: status character(len=*), intent(inout) :: gol ! ! !REVISION HISTORY: ! 12 Aug 2008 - Michael Kahnert, SMHI ! 6 Feb 2011 - Achim Strunk - ! ! !REMARKS: ! - see remarks of OPTICS_GET subroutine about PURE routine ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/get_refr_idx' ! refractive indices complex :: m_SO4, m_BC, m_OC, m_SOA, m_SS, m_DU, m_water ! volume fractions real :: v_SO4, v_BC, v_OC, v_SOA, v_SS, v_DU, v_water, water_iv integer :: i_test real :: vtot, v2 complex :: m00, m0, m1, m2 real :: rpls, ipls ! --- begin -------------------------------- ! Get the refractive indices from the lookup-tables and put them into complex numbers. m_so4 = cmplx( wdep%n(1),wdep%k(1) ) ! H2-SO4 + NH4NO3 m_bc = cmplx( wdep%n(2),wdep%k(2) ) ! BC m_oc = cmplx( wdep%n(3),wdep%k(3) ) ! POM m_soa = cmplx( wdep%n(4),wdep%k(4) ) ! SOA m_ss = cmplx( wdep%n(5),wdep%k(5) ) ! SS m_du = cmplx( wdep%n(6),wdep%k(6) ) ! DU m_water = cmplx( wdep%n(7),wdep%k(7) ) ! Water status = 0 ! no mixing for mode=6,7: if(mode.ge.6)then m_eff=m_DU return endif ! compute volume fractions: v_SO4=0. v_BC=0. v_OC=0. v_SOA=0. v_SS=0. v_DU=0. v_water=0. vtot=0. ! Added sanity check (15-7-2010 - P. Le Sager) : Avoid negative water ! mixing ratio! ! The bruggeman logically assumes that v_water is between 0 and 1, but ! this is never checked in the call chain : ! ECEarth_Optics_Step -> calculate_aop -> get_refr_idx [here] -> ! Bruggeman ! We do it here, with a warning since it reflects a problem upstream: if(water.lt.0.0)then !write (gol,*)" WARNING - [Get_refr_idx] : negative relative humidity..." ; call goPr !write (gol,*)" WARNING - [Get_refr_idx] : .....set to 0" ; call goPr water_iv=0.0 else water_iv=water endif if(mode.le.4)then v_SO4=SO4/so4_density v_water=water_iv/rol vtot=vtot+v_SO4+v_water endif if(mode.le.5)then !v_OC=OC/pom_density v_SOA=SOA/soa_density vtot=vtot+v_SOA end if if(mode.ge.2.and.mode.le.5)then v_BC=BC/carbon_density v_OC=OC/pom_density vtot=vtot+v_BC+v_OC endif if(mode.ge.3.and.mode.le.4)then v_SS=SS/ss_density vtot=vtot+v_SS endif if(mode.ge.3.and.mode.ne.5)then v_DU=DU/dust_density vtot=vtot+v_DU endif ! If vtot is zero, we will get 0.0/0.0's causing NaNs. In that case, the ! refractive index does not matter and will be set to (1.0,1.0e-9). The ! reason not to take (1.0,0.0) is that someone with humour might take ! the logarithm of the imaginary part. Dust particles get their usual ! refractive index, because they already returned m_DU. But that does ! not matter, because there are zero aerosols in this case. If (vtot .le. 1e-20) Then m_eff = Cmplx(1.0,1.0e-9) Else v_SO4=v_SO4/vtot v_BC=v_BC/vtot v_OC=v_OC/vtot v_SOA=v_SOA/vtot v_SS=v_SS/vtot v_DU=v_DU/vtot v_water=v_water/vtot !----------------------------------------------------------------------- ! effective medium computations for each mode !----------------------------------------------------------------------- if(mode.eq.1)then ! Bruggeman mixing rule for SO4 OC and water: m1=m_SO4 m2=m_SOA vtot=v_SO4+v_SOA v2=v_SOA/vtot call Bruggeman(m1,m2,v2,m0) m1=m0 m2=m_water v2=v_water call Bruggeman(m1,m2,v2,m0) elseif(mode.eq.2)then ! iterative Bruggeman mixing rule for SO4, OC, and water: m1=m_SO4 m2=m_OC vtot=v_SO4+v_OC v2=v_OC/vtot call Bruggeman(m1,m2,v2,m00) m1=m00 m2=m_SOA vtot=vtot+v_SOA v2=v_SOA/vtot call Bruggeman(m1,m2,v2,m00) m1=m00 m2=m_water vtot=vtot+v_water v2=v_water/vtot call Bruggeman(m1,m2,v2,m00) ! Maxwell-Garnett mixing rule for BC inclusions: m1=m00 m2=m_BC v2=v_BC call Maxwell_Garnett(m1,m2,v2,m0) elseif(mode.eq.3.or.mode.eq.4)then ! iterative Bruggeman mixing rule for SO4, OC, SS, and water: m1=m_SO4 m2=m_OC vtot=v_SO4+v_OC if ( vtot < TINY( vtot ) ) then v2=0.0 else v2=v_OC/vtot end if call Bruggeman(m1,m2,v2,m00) m1=m00 m2=m_SOA vtot=vtot+v_SOA if ( vtot < TINY( vtot ) ) then v2=0.0 else v2=v_SOA/vtot end if call Bruggeman(m1,m2,v2,m00) m1=m00 m2=m_SS vtot=vtot+v_SS if ( vtot < TINY( vtot ) ) then v2=0.0 else v2=v_SS/vtot end if call Bruggeman(m1,m2,v2,m00) m1=m00 m2=m_water vtot=vtot+v_water v2=v_water/vtot call Bruggeman(m1,m2,v2,m00) ! iterative Maxwell-Garnett mixing rule for BC and dust ! inclusions: m1=m00 m2=m_BC vtot=vtot+v_BC if ( vtot < TINY( vtot ) ) then v2=0.0 else v2=v_BC/vtot end if call Maxwell_Garnett(m1,m2,v2,m00) m1=m00 m2=m_DU v2=v_DU call Maxwell_Garnett(m1,m2,v2,m0) elseif(mode.eq.5)then m1=m_SOA m2=m_OC vtot=v_SOA+v_OC v2=v_OC/vtot call Bruggeman(m1,m2,v2,m00) ! Maxwell-Garnett mixing rule for BC inclusions: m1=m00 m2=m_BC v2=v_BC call Maxwell_Garnett(m1,m2,v2,m0) endif m_eff = m0 End If ! Debug : trap for a NAN (13-7-2010 - P. Le Sager) rpls=real(m_eff) ipls=imag(m_eff) IF ((rpls.NE.rpls).or.(ipls.NE.ipls)) then status = 1 write (gol,'(" GET_REFR_IDX-NAN: ", 3(E16.4,2x),i4,2x,7(E16.4,2x))') rpls, ipls, vtot, mode,& & SO4,BC,OC,SOA,SS,DU,water endif END SUBROUTINE GET_REFR_IDX !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: BRUGGEMAN ! ! !DESCRIPTION: Compute effective refractive index of a mixture of 2 components ! by use of the Bruggeman mixing rule !\\ !\\ ! !INTERFACE: ! PURE SUBROUTINE BRUGGEMAN(m1,m2,v2,m0) ! ! !INPUT PARAMETERS: ! complex, intent(in) :: m1,m2 real, intent(in) :: v2 ! ! !OUTPUT PARAMETERS: ! complex, intent(out) :: m0 ! ! !REVISION HISTORY: ! 12 Aug 2008 - Michael Kahnert, SMHI ! 6 Feb 2011 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC complex m1s,m2s,mt real fac1,fac2 if(v2.eq.1.0)then m0=m2 return elseif(v2.eq.0.0)then m0=m1 return endif fac1=2.-3.*v2 fac2=3.*v2-1. m1s=m1**2 m2s=m2**2 mt=m1s*fac1+m2s*fac2 m0=1./16.*mt**2+0.5*m1s*m2s m0=csqrt(m0) m0=m0+0.25*mt m0=csqrt(m0) END SUBROUTINE BRUGGEMAN !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: MAXWELL_GARNETT ! ! !DESCRIPTION: Compute effective refractive index for a medium consisting of ! a matrix with refractive index m1 and inclusions with refractive ! index m2 and volume fraction v2 by use of the Maxwell-Garnett ! mixing rule. !\\ !\\ ! !INTERFACE: ! PURE SUBROUTINE MAXWELL_GARNETT( m1, m2, v2, m0) ! ! !INPUT PARAMETERS: ! complex, intent(in) :: m1, m2 real, intent(in) :: v2 ! ! !OUTPUT PARAMETERS: ! complex, intent(out) :: m0 ! ! !REVISION HISTORY: ! 12 Aug 2008 - Michael Kahnert, SMHI ! 6 Feb 2011 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC complex :: m1s,m2s real :: fac1,fac2 fac1=3.0-2.0*v2 fac2=3.0-v2 m1s=m1**2 m2s=m2**2 m0=m2s*(fac1*m1s+2.0*v2*m2s)/(v2*m1s+fac2*m2s) m0=csqrt(m0) END SUBROUTINE MAXWELL_GARNETT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OPTICS_CALCULATE_AOP ! ! !DESCRIPTION: This routine writes on aop_out_* (module wide parameters) ! the retrieved aerosol properties. The caller has to ensure ! that these fields have been allocated properly. ! IMPORTANT: OC is actually POM. ! Remember converting OC to POM when sending it to this method. !\\ !\\ ! !INTERFACE: ! SUBROUTINE OPTICS_CALCULATE_AOP( nwl, wdep, lvec, ecearth_units, & aop_out_ext, aop_out_a, aop_out_g, status, aop_out_add ) ! ! !USES: ! use binas, only : rol, twopi use chem_param, only : ss_density, dust_density, carbon_density use chem_param, only : pom_density, so4_density, soa_density Use mo_aero_m7, only : nmod, cmedr2mmedr, sigma ! ! !INPUT PARAMETERS: ! integer, intent(in) :: nwl type(wavelendep), intent(in), dimension(nwl) :: wdep ! wavelength depending properties integer, intent(in) :: lvec logical, intent(in) :: ecearth_units ! ! !OUTPUT PARAMETERS: ! Real, Dimension(:,:,:), intent(out) :: aop_out_ext ! extinctions Real, Dimension(:,:), intent(out) :: aop_out_a ! single scattering albedo Real, Dimension(:,:), intent(out) :: aop_out_g ! assymetry factor integer, intent(out) :: status Real, Dimension(:,:,:), intent(out), optional :: aop_out_add ! additional parameters ! ! !REVISION HISTORY: ! 12 Aug 2008 - Michael Kahnert, SMHI ! 6 Feb 2011 - Achim Strunk - ! 23 Jan 2013 - Narcisa Banda - included OMP (PLS: commented in TM5-MP) ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC real, dimension(:), allocatable :: NCsca, incext real :: Cexti, ai, gi, NCscai, xg complex :: m_eff real, dimension(:), allocatable :: lnsigma integer :: i, imode, ivec!, lss, lee integer :: statusomp Logical :: coarse real :: totvoldry, modfrac ! real :: t1, t2 , omp_get_wtime Real, Dimension(:,:,:), Pointer :: cext_table, a_table, g_table ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/optics_calculate_AOP' ! --- begin -------------------------------- allocate( NCsca ( lvec ) ) ! allocate( NCscai( lvec ) ) ! allocate( Cexti ( lvec ) ) ! allocate( ai ( lvec ) ) ! allocate( gi ( lvec ) ) ! allocate( xg ( lvec ) ) ! allocate( m_eff ( lvec ) ) allocate( incext( lvec ) ) !t1= omp_get_wtime() status = 0 !pls! !$OMP PARALLEL & !pls! !$OMP default (none) & !pls! !$OMP shared (nwl, lvec, NCsca, incext, & !pls! !$OMP wdep, aop_in, aop_out_Ext, aop_out_g, aop_out_a, aop_out_add, sigma, cmedr2mmedr,& !pls! !$OMP cext_200, a_200, g_200, cext_159, a_159, g_159, status, gol) & !pls! !$OMP private (lnsigma, i, imode, modfrac, coarse, totvoldry, lss, lee, ivec, & !pls! !$OMP NCscai, Cexti, ai, gi, xg, m_eff, cext_table, a_table, g_table, statusomp ) !pls! call my_omp_domain (1, lvec, lss, lee) ! Sulphate based on OPAC (Hess et al., 1998): !======================================================================= ! Get refractive indices of each component at the given wavelength: !======================================================================= do i = 1, nwl ! loop over wavelengths if( wdep(i)%split .or. wdep(i)%insitu ) then allocate( lnsigma( nmod ) ) lnsigma = log(sigma) end if aop_out_Ext(:,i,:) = 0.0 aop_out_g (:,i) = 0.0 aop_out_a (:,i) = 0.0 if( wdep(i)%insitu ) aop_out_add(:,i,:) = 0.0 NCsca (:) = 0.0 do imode = 1, 7 ! loop over M7 modes coarse = (imode .eq. 4 .or. imode .eq. 7) If (coarse) then cext_table => cext_200 a_table => a_200 g_table => g_200 Else cext_table => cext_159 a_table => a_159 g_table => g_159 End if !======================================================================= ! Compute effective refractive index of internally mixed aerosols ! for each grid cell and mode: !======================================================================= do ivec = 1, lvec !write(1001,*) aop_in(ivec)%SO4(imode), aop_in(ivec)%NO3(imode) !write(1002,*) aop_in(ivec)%SOA(imode) !write(1003,*) aop_in(ivec)%BC(imode), aop_in(ivec)%OC(imode) !write(1004,*) aop_in(ivec)%SS(imode), aop_in(ivec)%DU(imode) call get_refr_idx( wdep(i), & aop_in(ivec)%SO4(imode) + aop_in(ivec)%NO3(imode), & ! H2-SO4 + NH4NO3 aop_in(ivec)%BC (imode), aop_in(ivec)%OC (imode), & aop_in(ivec)%SOA (imode), & aop_in(ivec)%SS (imode), aop_in(ivec)%DU (imode), & aop_in(ivec)%h2o(imode), imode, m_eff, statusomp, gol) if (statusomp ==1) then ! Problem tracked to emilnox, should not happen anymore ! but keep the printout and hack for now !FIX to finish crescendo write (gol,'(" Problem with GET_REFR_IDX-NAN inputs: ", 9(E16.4,2x))') & aop_in(ivec)%SO4(imode), aop_in(ivec)%NO3(imode), aop_in(ivec)%h2o(imode), & aop_in(ivec)%SOA(imode),aop_in(ivec)%BC(imode), aop_in(ivec)%OC(imode),& aop_in(ivec)%SS(imode), aop_in(ivec)%DU(imode) call GoErr ! Assume edge case not caught. Set m_eff to default (1.0,1.0e-9), same as in ref_idx ! (subroutine problems with inputs SO4 and WATER?) status = 0 m_eff= Cmplx(1.0,1.0e-9) ! If all edge cases were handled correctly, that should be: !status = 1 !IF_ERROR_RETURN( status=1 ) endif ! cmk added towpi for new netcdf lookup table xg = twopi*aop_in(ivec)%rg(imode) / wdep(i)%wl !======================================================================= ! get aerosol optical properties from data base for each mode !======================================================================= ! add extra safeguard for negative xg. It is unphysical, but have a safeguard anyway ! added for refr_idx_nan problem if (xg .gt. 0.0) then call Optics_Get(m_eff, xg, Cexti, ai, gi, cext_table, a_table, g_table, statusomp, gol ) else Cexti=0.0 ai=1.0 gi=1.0 end if if (statusomp ==1) then call GoErr status = 1 IF_ERROR_RETURN( status=1 ) endif ! Multiply Cext with lambda^2 to get the cross section. Cexti = Cexti*(wdep(i)%wl**2) ! this here is extinction coefficient in this mode incext(ivec) = aop_in(ivec)%numdens(imode) * Cexti ! sum up partial coefficients Aop_out_ext(ivec,i,1) = Aop_out_Ext(ivec,i,1) + incext(ivec) ! scattering portion NCscai = ai * incext(ivec) ! sum up weights for average (both albedo and asymmetry) NCsca (ivec) = NCsca (ivec) + NCscai aop_out_g(ivec,i) = aop_out_g(ivec,i) + NCscai * gi end do ! Split extinction to separate contributions from constituents in modes. ! A volume mixing is assumed (in contrast to the explicit mixing in get_refr_ind). if( wdep(i)%split .or. wdep(i)%insitu) then do ivec = 1, lvec if (wdep(i)%split) then ! The fine-mode contributions to the extinction ! includes the contributions from particles ! with (wet) diameters smaller than 1 micron. ! For wet particles, only part of the accumulation mode ! should be included, and the coarse mode should be excluded. if (.not.coarse) then ! Currently, the contribution of the accumulation mode ! is approximated using weighting by volume scaling factor modfrac. ! For extinction, area weighting would probably be more appropriate!! if (imode .eq. 3 .or. imode .eq. 6 ) then ! - convert number mean radius to volume mean radius (by cmedr2mmedr(imode)) ! - 1 micron diameter --> radius is 0.5 microns (rg is also in microns) modfrac = 0.5 + 0.5 * ERF( log( 0.5 / (aop_in(ivec)%rg(imode) * cmedr2mmedr(imode)) ) / & (sqrt(2.0) * lnsigma(imode)) ) else ! Include full nucleation and Aitken mode contributions. modfrac = 1.0 endif aop_out_ext(ivec,i,10) = aop_out_ext(ivec,i,10) + modfrac * incext(ivec) endif ! total volume from so4/no3/bc/oc/ss/du in this mode (ATTENTION: DRY!!) ! take no3 as so4 totvoldry = aop_in(ivec)%so4(imode)/so4_density + aop_in(ivec)%no3(imode)/so4_density + & aop_in(ivec)%bc (imode)/carbon_density + aop_in(ivec)%oc (imode)/pom_density + & aop_in(ivec)%soa (imode)/soa_density + & aop_in(ivec)%ss (imode)/ss_density + aop_in(ivec)%du (imode)/dust_density ! check whether there is some volume available ! otherwise assign zeros to extinction increments if( totvoldry < tiny(totvoldry) ) then write(gol,'("WARNING: no volume in mode (",i3,"). assigning zero extinctions")') imode; call goPr cycle end if ! H2-SO4 contribution aop_out_ext(ivec,i,2) = aop_out_ext(ivec,i,2) + incext(ivec) * (aop_in(ivec)%so4(imode)/so4_density ) / totvoldry ! BC contribution aop_out_ext(ivec,i,3) = aop_out_ext(ivec,i,3) + incext(ivec) * (aop_in(ivec)%bc (imode)/carbon_density) / totvoldry ! POM contribution aop_out_ext(ivec,i,4) = aop_out_ext(ivec,i,4) + incext(ivec) * (aop_in(ivec)%oc (imode)/pom_density ) / totvoldry ! SOA contribution aop_out_ext(ivec,i,5) = aop_out_ext(ivec,i,5) + incext(ivec) * (aop_in(ivec)%soa (imode)/soa_density ) / totvoldry ! SS contribution aop_out_ext(ivec,i,6) = aop_out_ext(ivec,i,6) + incext(ivec) * (aop_in(ivec)%ss (imode)/ss_density ) / totvoldry ! DU contribution aop_out_ext(ivec,i,7) = aop_out_ext(ivec,i,7) + incext(ivec) * (aop_in(ivec)%du (imode)/dust_density ) / totvoldry ! NH4NO3 contribution aop_out_ext(ivec,i,8) = aop_out_ext(ivec,i,8) + incext(ivec) * (aop_in(ivec)%no3(imode)/so4_density ) / totvoldry ! Fine-mode contributions for dust and sea salt if (.not.coarse) then aop_out_ext(ivec,i,11) = aop_out_ext(ivec,i,11) + incext(ivec) * (aop_in(ivec)%du (imode)/dust_density ) / totvoldry * modfrac aop_out_ext(ivec,i,12) = aop_out_ext(ivec,i,12) + incext(ivec) * (aop_in(ivec)%ss (imode)/ss_density ) / totvoldry * modfrac endif endif ! Water contribution: ! Get optical properties without water, the difference will be extinction due to ! water existence ! - mis-use gi for this gi = 0.0 call get_refr_idx( wdep(i), & aop_in(ivec)%SO4(imode) + aop_in(ivec)%NO3(imode), & aop_in(ivec)%BC (imode), aop_in(ivec)%OC (imode), & aop_in(ivec)%SOA (imode), & aop_in(ivec)%SS (imode), aop_in(ivec)%DU (imode), & gi, imode, m_eff, statusomp, gol) if (statusomp ==1) then write (gol,'("Problem with GET_REFR_IDX-NAN inputs: ", 9(E16.4,2x))') & aop_in(ivec)%SO4(imode), aop_in(ivec)%NO3(imode), aop_in(ivec)%h2o(imode), & aop_in(ivec)%SOA(imode),aop_in(ivec)%BC(imode), aop_in(ivec)%OC(imode),& aop_in(ivec)%SS(imode), aop_in(ivec)%DU(imode) call GoErr ! Assume edge case not caught. Set m_eff to default (1.0,1.0e-9), same as in ref_idx ! (subroutine problems with inputs SO4 and WATER?) status = 0 m_eff= Cmplx(1.0,1.0e-9) ! If all edge cases were handled correctly, that should be: !status = 1 !IF_ERROR_RETURN( status=1 ) !status = 1 endif ! here we need the dry radius!!! !>>> TvN ! 2*pi should be included (see comment above) !xg = aop_in(ivec)%rgd(imode) / wdep(i)%wl xg = twopi*aop_in(ivec)%rgd(imode) / wdep(i)%wl !<<< TvN ! TB ! add extra safeguard for negative xg, unphysical but have a safeguard ! added for refr_idx_nan problem if (xg .gt. 0.0) then call Optics_Get(m_eff, xg, Cexti, ai, gi, cext_table, a_table, g_table, statusomp, gol ) else Cexti=0.0 ai=1.0 gi=1.0 end if if (statusomp ==1) then call GoErr status = 1 IF_ERROR_RETURN( status=1 ) endif Cexti = Cexti*(wdep(i)%wl**2) if (wdep(i)%split) then ! add difference to water subarray Aop_out_ext(ivec,i,9) = aop_out_ext(ivec,i,9) + (incext(ivec) - aop_in(ivec)%numdens(imode) * Cexti) endif if (wdep(i)%insitu) then ! Surface dry extinction and absorption: !>>> TvN ! Remove cut off for the total values: !modfrac = 0.5 + 0.5 * ERF( log( 5.0 / (aop_in(ivec)%rgd(imode) * cmedr2mmedr(imode)) ) / & ! (sqrt(2.0) * lnsigma(imode)) ) ! extinction and absorption (extinction-scattering): !Aop_out_add(ivec,i,1) = aop_out_add(ivec,i,1) + modfrac * aop_in(ivec)%numdens(imode) * Cexti !Aop_out_add(ivec,i,2) = aop_out_add(ivec,i,2) + modfrac * aop_in(ivec)%numdens(imode) * Cexti * (1. - ai) aop_out_add(ivec,i,1) = aop_out_add(ivec,i,1) + aop_in(ivec)%numdens(imode) * Cexti aop_out_add(ivec,i,2) = aop_out_add(ivec,i,2) + aop_in(ivec)%numdens(imode) * Cexti * (1. - ai) aop_out_add(ivec,i,3) = aop_out_add(ivec,i,3) + aop_in(ivec)%numdens(imode) * Cexti * ai * gi ! Add fine-mode contributions ! For dry aerosol, the fine-mode optical properties include ! the full accumulation mode, but not the coarse mode: if (.not.coarse) then aop_out_add(ivec,i,4) = aop_out_add(ivec,i,4) + aop_in(ivec)%numdens(imode) * Cexti aop_out_add(ivec,i,5) = aop_out_add(ivec,i,5) + aop_in(ivec)%numdens(imode) * Cexti * (1. - ai) endif !<<< TvN endif end do ! ivec end if enddo ! modes if (ecearth_units) then ! return extinction due to absorption ! and asymmetry factor multiplied by extinction due to scattering ! and convert um**2/m**3 into 1/m (as for extinction below) Aop_out_a(:,i) = ( Aop_out_Ext(:,i,1) - NCsca(:) ) * 1.e-12 Aop_out_g(:,i) = Aop_out_g(:,i) * 1.e-12 else ! return single-scattering albedo and asymmetry factor ! take into account small extinction values where( Aop_out_Ext(:,i,1) > tiny(0.0) ) Aop_out_a(:,i) = NCsca(:) / Aop_out_Ext(:,i,1) elsewhere ! No extinction -> Fill 1.0 for albedo because it looks clean. Aop_out_a(:,i) = 1.0 endwhere ! take into account small extinction values where( NCsca(:) > tiny(0.0) ) Aop_out_g(:,i) = Aop_out_g(:,i) / NCsca(:) elsewhere ! No scattering at all -> Fill 1.0 for asymmetry because it looks clean. Aop_out_g(:,i) = 1.0 endwhere endif ! convert um**2/m**3 into 1/m Aop_out_Ext(:,i,:) = Aop_out_Ext(:,i,:) * 1e-12 if( wdep(i)%insitu) then ! take into account small extinction values where( aop_out_add(:,i,1) - aop_out_add(:,i,2) > tiny(0.0) ) aop_out_add(:,i,3) = aop_out_add(:,i,3) / (aop_out_add(:,i,1)-aop_out_add(:,i,2)) elsewhere ! No scattering at all -> Fill 1.0 for asymmetry because it looks clean. aop_out_add(:,i,3) = 1.0 endwhere Aop_out_add(:,i,1) = aop_out_add(:,i,1) * 1e-12 Aop_out_add(:,i,2) = aop_out_add(:,i,2) * 1e-12 Aop_out_add(:,i,4) = aop_out_add(:,i,4) * 1e-12 Aop_out_add(:,i,5) = aop_out_add(:,i,5) * 1e-12 endif if( wdep(i)%split .or. wdep(i)%insitu ) then deallocate( lnsigma ) endif !======================================================================= enddo ! loop over wavelengths Nullify(Cext_table) Nullify(a_table) Nullify(g_table) !$OMP END PARALLEL if (status==1) then call goPr IF_NOTOK_RETURN(status=1) endif !do imode = 1,7 ! print *, 'Radius,mode :', imode, sum(aop_in(:)%rg(imode))/lvec ! print *, 'numden,mode :', imode, sum(aop_in(:)%numdens(imode))/lvec !enddo deallocate( NCsca, incext ) !do i = 1,nwl ! print *, 'AOD per grid box:', wdep(i)%wl, sum(Aop_out_Ext(:,i,1))/lvec !enddo ! ok status = 0 END SUBROUTINE OPTICS_CALCULATE_AOP !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OPTICS_AOP_GET ! ! !DESCRIPTION: Initialise the fields "aop_in" and then calculate the ! optical properties through a call to optics_calculate_aop. !\\ !\\ ! !INTERFACE: ! SUBROUTINE OPTICS_AOP_GET( lvec, region, nwav, wdep, ncontr, ecearth_units, & aop_out_ext, aop_out_a, aop_out_g, status, aop_out_add ) ! ! !USES: ! USE TM5_DISTGRID, ONLY : dgrid, Get_DistGrid use Meteo, only : Set Use Meteodata, only : gph_dat, m_dat Use Dims, only : im, jm, lm use global_data, only : region_dat, mass_dat Use chem_param, only : inus_n, iso4nus, isoanus, nmod use chem_param, only : iais_n, iso4ais, ibcais, ipomais, isoaais Use chem_param, only : iacs_n, iso4acs, ibcacs, ipomacs, isoaacs, issacs, iduacs, ino3_a use chem_param, only : imsa use chem_param, only : icos_n, iso4cos, ibccos, ipomcos, isoacos, isscos, iducos use chem_param, only : iaii_n, ibcaii, ipomaii, isoaaii use chem_param, only : iaci_n, icoi_n, iduaci, iducoi use chem_param, only : h2so4_factor, nh4no3_factor use chem_param, only : so4_density, nh4no3_density, msa_density Use m7_data, only : h2o_mode, rw_mode, rwd_mode ! ! !INPUT PARAMETERS: ! Integer, Intent(In) :: region, lvec Integer, Intent(In) :: nwav, ncontr type(wavelendep), dimension(nwav), Intent(In) :: wdep logical, intent(in) :: ecearth_units ! ! !OUTPUT PARAMETERS: ! real, dimension(lvec,nwav,ncontr), intent(out) :: aop_out_ext ! extinctions real, dimension(lvec,nwav), intent(out) :: aop_out_a ! single scattering albedo real, dimension(lvec,nwav), intent(out) :: aop_out_g ! assymetry factor integer, intent(out) :: status real, dimension(:,:,:), intent(out), optional :: aop_out_add ! additional parameters ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - v0 ! 24 Jun 2011 - Achim Strunk - Added NO3 to allow for explicit ! AOD split of (SO4+NO3) ! 20 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/optics_aop_get' real, dimension(:,:,:,:), pointer :: rm real, dimension(:,:,:), pointer :: m real, dimension(:,:,:), allocatable :: volume Integer :: i, j, l, iloop, i1, i2, j1, j2, lmr ! --- start ------------------------------ ! ensure met fields being available call Set( gph_dat(region), status, used=.true. ) ! air & tracers mass m => m_dat(region)%data rm => mass_dat(region)%rm ! tile grid size CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! Generalize to allow for reduced number of levels ! when called from ECEarth_Optics_Step !lmr=lm(region) lmr = lvec / ( (i2-i1+1)*(j2-j1+1) ) allocate( aop_in(lvec) ) do i = 1, nmod aop_in(:)%so4 (i) = 0.0 ; aop_in(:)%bc (i) = 0.0 aop_in(:)%oc (i) = 0.0 ; aop_in(:)%soa (i) = 0.0 aop_in(:)%ss (i) = 0.0 aop_in(:)%du (i) = 0.0 ; aop_in(:)%h2o (i) = 0.0 aop_in(:)%numdens(i) = 0.0 ; aop_in(:)%rg (i) = 0.0 aop_in(:)%rgd (i) = 0.0 ; aop_in(:)%no3 (i) = 0.0 end do !>>> TvN ! In M7 sulphate is assumed to be H2-SO4 with corresponding particle density so4_density ! The sulphate mass should therefore also include the small contribution from the H atoms ! NUS aop_in(:)%so4(1) = reshape( 1.e9 * h2so4_factor * rm(i1:i2,j1:j2,1:lmr,iso4nus) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%soa(1) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoanus) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! AIS aop_in(:)%so4(2) = reshape( 1.e9 * h2so4_factor * rm(i1:i2,j1:j2,1:lmr,iso4ais) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%bc (2) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ibcais ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%oc (2) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ipomais) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%soa(2) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoaais) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! ACS (additional: NO3) ! The contribution from methane sulfonate (MSA-) aerosol is added ! to that for sulfate. ! As the addition is done by volume, ! we need to account for the difference in densities ! (as done below for ammonium nitrate). aop_in(:)%so4(3) = reshape( 1.e9 * ( h2so4_factor * rm(i1:i2,j1:j2,1:lmr,iso4acs) + & (so4_density / msa_density) * rm(i1:i2,j1:j2,1:lmr,imsa) ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! Since nh4no3_density is the density of NH4NO3, the contribution from NH4 should be included. ! Moreover, assuming the same refractive index for NH4NO3 as for H2-SO4, ! the contributions from both components can be added by volume; ! thus we need to account for the difference in densities. ! Estimates of the refractive index of NH4NO3 are available from literature ! (e.g. Lowenthal et al., Atmos. Environ., 2000). ! For practical purposes, it can be set equal to the value used for sulfate, ! i.e. the value for a solution containing 75% H2SO4 (Fenn et al., 1985). aop_in(:)%no3(3) = reshape( 1.e9 * nh4no3_factor * (so4_density / nh4no3_density) * & rm(i1:i2,j1:j2,1:lmr,ino3_a ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%bc (3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ibcacs ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%oc (3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ipomacs) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%soa(3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoaacs) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%ss (3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,issacs ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%du (3) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,iduacs ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! COS aop_in(:)%so4(4) = reshape( 1.e9 * h2so4_factor * rm(i1:i2,j1:j2,1:lmr,iso4cos) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) !<<< TvN aop_in(:)%bc (4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ibccos ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%oc (4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ipomcos) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%soa(4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoacos) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%ss (4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isscos ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%du (4) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,iducos ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! AII aop_in(:)%bc (5) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ibcaii ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%oc (5) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,ipomaii) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%soa(5) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,isoaaii) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! ACI aop_in(:)%du (6) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,iduaci ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! COI aop_in(:)%du (7) = reshape( 1.e9 * rm(i1:i2,j1:j2,1:lmr,iducoi ) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! Water in (hydrophillic) modes aop_in(:)%h2o(1) = reshape( 1.e9 * h2o_mode(region,1)%d3(i1:i2,j1:j2,1:lmr) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%h2o(2) = reshape( 1.e9 * h2o_mode(region,2)%d3(i1:i2,j1:j2,1:lmr) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%h2o(3) = reshape( 1.e9 * h2o_mode(region,3)%d3(i1:i2,j1:j2,1:lmr) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%h2o(4) = reshape( 1.e9 * h2o_mode(region,4)%d3(i1:i2,j1:j2,1:lmr) / & m(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rg (1) = reshape( 1.e6 * rw_mode (region,1)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rg (2) = reshape( 1.e6 * rw_mode (region,2)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rg (3) = reshape( 1.e6 * rw_mode (region,3)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rg (4) = reshape( 1.e6 * rw_mode (region,4)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rg (5) = reshape( 1.e6 * rw_mode (region,5)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rg (6) = reshape( 1.e6 * rw_mode (region,6)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rg (7) = reshape( 1.e6 * rw_mode (region,7)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) ! dry radius for soluble modes / rest equals the usual radii aop_in(:)%rgd(1) = reshape( 1.e6 * rwd_mode(region,1)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rgd(2) = reshape( 1.e6 * rwd_mode(region,2)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rgd(3) = reshape( 1.e6 * rwd_mode(region,3)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rgd(4) = reshape( 1.e6 * rwd_mode(region,4)%d3(i1:i2,j1:j2,1:lmr), (/lvec/) ) aop_in(:)%rgd(5) = aop_in(:)%rg (5) aop_in(:)%rgd(6) = aop_in(:)%rg (6) aop_in(:)%rgd(7) = aop_in(:)%rg (7) ! volume for number density allocate( volume( i1:i2, j1:j2, 1:lmr ) ) do j = j1, j2 volume(:,j,:) = ( gph_dat(region)%data(i1:i2, j, 2:lmr+1) - & gph_dat(region)%data(i1:i2, j, 1:lmr ) ) * & region_dat(region)%dxyp(j) ! m3 end do aop_in(:)%numdens(1) = reshape( rm(i1:i2,j1:j2,1:lmr,inus_n) / volume, (/lvec/) ) aop_in(:)%numdens(2) = reshape( rm(i1:i2,j1:j2,1:lmr,iais_n) / volume, (/lvec/) ) aop_in(:)%numdens(3) = reshape( rm(i1:i2,j1:j2,1:lmr,iacs_n) / volume, (/lvec/) ) aop_in(:)%numdens(4) = reshape( rm(i1:i2,j1:j2,1:lmr,icos_n) / volume, (/lvec/) ) aop_in(:)%numdens(5) = reshape( rm(i1:i2,j1:j2,1:lmr,iaii_n) / volume, (/lvec/) ) aop_in(:)%numdens(6) = reshape( rm(i1:i2,j1:j2,1:lmr,iaci_n) / volume, (/lvec/) ) aop_in(:)%numdens(7) = reshape( rm(i1:i2,j1:j2,1:lmr,icoi_n) / volume, (/lvec/) ) deallocate( volume ) ! check valid ranges in particle sizes (might be zero) where( aop_in%rg (1) .lt. 1.E-15 ) aop_in%rg (1) = 1.E-15 where( aop_in%rgd(1) .lt. 1.E-15 ) aop_in%rgd(1) = 1.E-15 where( aop_in%rg (2) .lt. 1.E-15 ) aop_in%rg (2) = 1.E-15 where( aop_in%rgd(2) .lt. 1.E-15 ) aop_in%rgd(2) = 1.E-15 where( aop_in%rg (3) .lt. 1.E-15 ) aop_in%rg (3) = 1.E-15 where( aop_in%rgd(3) .lt. 1.E-15 ) aop_in%rgd(3) = 1.E-15 where( aop_in%rg (4) .lt. 1.E-15 ) aop_in%rg (4) = 1.E-15 where( aop_in%rgd(4) .lt. 1.E-15 ) aop_in%rgd(4) = 1.E-15 where( aop_in%rg (5) .lt. 1.E-15 ) aop_in%rg (5) = 1.E-15 where( aop_in%rgd(5) .lt. 1.E-15 ) aop_in%rgd(5) = 1.E-15 where( aop_in%rg (6) .lt. 1.E-15 ) aop_in%rg (6) = 1.E-15 where( aop_in%rgd(6) .lt. 1.E-15 ) aop_in%rgd(6) = 1.E-15 where( aop_in%rg (7) .lt. 1.E-15 ) aop_in%rg (7) = 1.E-15 where( aop_in%rgd(7) .lt. 1.E-15 ) aop_in%rgd(7) = 1.E-15 nullify(rm) nullify(m) ! aop_out_ext=0.0 aop_out_a=0.0 aop_out_g=0.0 if (present(aop_out_add)) then call optics_calculate_aop( nwav, wdep, lvec, ecearth_units, & aop_out_ext, aop_out_a, aop_out_g, status, aop_out_add ) else call optics_calculate_aop( nwav, wdep, lvec, ecearth_units, & aop_out_ext, aop_out_a, aop_out_g, status ) endif IF_NOTOK_RETURN(status=1) ! OK Deallocate( aop_in ) status = 0 END SUBROUTINE OPTICS_AOP_GET !EOC END MODULE OPTICS