Index: LMDZ5/trunk/DefLists/field_def_lmdz.xml =================================================================== --- LMDZ5/trunk/DefLists/field_def_lmdz.xml (revision 2689) +++ LMDZ5/trunk/DefLists/field_def_lmdz.xml (revision 2690) @@ -613,4 +613,29 @@ + + + + + + + + + + + + + + + + + + + + + + + + + Index: LMDZ5/trunk/DefLists/file_def_histStratAer_lmdz.xml =================================================================== --- LMDZ5/trunk/DefLists/file_def_histStratAer_lmdz.xml (revision 2690) +++ LMDZ5/trunk/DefLists/file_def_histStratAer_lmdz.xml (revision 2690) @@ -0,0 +1,35 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Index: LMDZ5/trunk/bld.cfg =================================================================== --- LMDZ5/trunk/bld.cfg (revision 2689) +++ LMDZ5/trunk/bld.cfg (revision 2690) @@ -29,4 +29,5 @@ src::rrtm %RRTM src::dust %DUST +src::strataer %STRATAER src::grid %SRC_PATH/grid src::filtrez %FILTRE @@ -108,2 +109,3 @@ bld::tool::SHELL /bin/bash bld::tool::SHELL /bin/ksh +bld::tool::SHELL /bin/ksh Index: LMDZ5/trunk/libf/dyn3d_common/infotrac.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3d_common/infotrac.F90 (revision 2689) +++ LMDZ5/trunk/libf/dyn3d_common/infotrac.F90 (revision 2690) @@ -41,18 +41,29 @@ CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym - ! CRisi: cas particulier des isotopes - LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso - INTEGER :: niso_possibles - PARAMETER ( niso_possibles=5) - real, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal - LOGICAL, DIMENSION(niso_possibles),SAVE :: use_iso - INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase) - INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot - INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot - INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numéro de la zone de tracage en fn de nqtot - INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numéro de la zone de tracage en fn de nqtot - INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles - INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: index_trac ! numéro ixt en fn izone, indnum entre 1 et niso - INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso +! CRisi: cas particulier des isotopes + LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso + INTEGER :: niso_possibles + PARAMETER ( niso_possibles=5) + REAL, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal + LOGICAL, DIMENSION(niso_possibles),SAVE :: use_iso + INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: iqiso ! donne indice iq en fn de (ixt,phase) + INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot + INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot + INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: zone_num ! donne numéro de la zone de tracage en fn de nqtot + INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: phase_num ! donne numéro de la zone de tracage en fn de nqtot + INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles + INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE :: index_trac ! numéro ixt en fn izone, indnum entre 1 et niso + INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso + +#ifdef CPP_StratAer +!--CK/OB for stratospheric aerosols + INTEGER, SAVE :: nbtr_bin + INTEGER, SAVE :: nbtr_sulgas + INTEGER, SAVE :: id_OCS_strat + INTEGER, SAVE :: id_SO2_strat + INTEGER, SAVE :: id_H2SO4_strat + INTEGER, SAVE :: id_BIN01_strat + INTEGER, SAVE :: id_TEST_strat +#endif CONTAINS @@ -141,4 +152,10 @@ CALL abort_gcm('infotrac_init','You must compile with cpp key REPROBUS',1) #endif + ELSE IF (type_trac == 'coag') THEN + WRITE(lunout,*) 'Tracers are treated for COAGULATION tests : type_trac=', type_trac +#ifndef CPP_StratAer + WRITE(lunout,*) 'To run this option you must add cpp key StratAer and compile with StratAer code' + CALL abort_gcm('infotrac_init','You must compile with cpp key StratAer',1) +#endif ELSE IF (type_trac == 'lmdz') THEN WRITE(lunout,*) 'Tracers are treated in LMDZ only : type_trac=', type_trac @@ -148,5 +165,4 @@ END IF - ! Test if config_inca is other then none for run without INCA IF (type_trac/='inca' .AND. config_inca/='none') THEN @@ -155,5 +171,4 @@ END IF - !----------------------------------------------------------------------- ! @@ -162,5 +177,5 @@ ! !----------------------------------------------------------------------- - IF (type_trac == 'lmdz' .OR. type_trac == 'repr') THEN + IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag') THEN OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr) IF(ierr.EQ.0) THEN @@ -171,10 +186,10 @@ WRITE(lunout,*) trim(modname),': Problem in opening traceur.def' WRITE(lunout,*) trim(modname),': WARNING using defaut values' - if (planet_type=='earth') then + IF (planet_type=='earth') THEN nqtrue=4 ! Default value for Earth - else + ELSE nqtrue=1 ! Default value for other planets - endif - END IF + ENDIF + ENDIF !jyg< !! if ( planet_type=='earth') then @@ -211,5 +226,5 @@ ALLOCATE(hadv_inca(nbtr), vadv_inca(nbtr)) - END IF ! type_trac + ENDIF ! type_trac !>jyg @@ -266,5 +281,5 @@ ! Get choice of advection schema from file tracer.def or from INCA !--------------------------------------------------------------------- - IF (type_trac == 'lmdz' .OR. type_trac == 'repr') THEN + IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag') THEN IF(ierr.EQ.0) THEN ! Continue to read tracer.def @@ -346,20 +361,58 @@ END DO - if ( planet_type=='earth') then + IF ( planet_type=='earth') THEN !CR: nombre de traceurs de l eau - if (tnom_0(3) == 'H2Oi') then + IF (tnom_0(3) == 'H2Oi') THEN nqo=3 - else + ELSE nqo=2 - endif + ENDIF ! For Earth, water vapour & liquid tracers are not in the physics nbtr=nqtrue-nqo - else + ELSE ! Other planets (for now); we have the same number of tracers ! in the dynamics than in the physics nbtr=nqtrue - endif - - ENDIF ! (type_trac == 'lmdz' .OR. type_trac == 'repr') + ENDIF + +#ifdef CPP_StratAer + IF (type_trac == 'coag') THEN + nbtr_bin=0 + nbtr_sulgas=0 + DO iq=1,nqtrue + IF (tnom_0(iq)(1:3)=='BIN') THEN !check if tracer name contains 'BIN' + nbtr_bin=nbtr_bin+1 + ENDIF + IF (tnom_0(iq)(1:3)=='GAS') THEN !check if tracer name contains 'GAS' + nbtr_sulgas=nbtr_sulgas+1 + ENDIF + ENDDO + print*,'nbtr_bin=',nbtr_bin + print*,'nbtr_sulgas=',nbtr_sulgas + DO iq=1,nqtrue + IF (tnom_0(iq)=='GASOCS') THEN + id_OCS_strat=iq-nqo + ENDIF + IF (tnom_0(iq)=='GASSO2') THEN + id_SO2_strat=iq-nqo + ENDIF + IF (tnom_0(iq)=='GASH2SO4') THEN + id_H2SO4_strat=iq-nqo + ENDIF + IF (tnom_0(iq)=='BIN01') THEN + id_BIN01_strat=iq-nqo + ENDIF + IF (tnom_0(iq)=='GASTEST') THEN + id_TEST_strat=iq-nqo + ENDIF + ENDDO + print*,'id_OCS_strat =',id_OCS_strat + print*,'id_SO2_strat =',id_SO2_strat + print*,'id_H2SO4_strat=',id_H2SO4_strat + print*,'id_BIN01_strat=',id_BIN01_strat + ENDIF +#endif + + ENDIF ! (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac = 'coag') !jyg< ! Index: LMDZ5/trunk/libf/phylmd/StratAer/aer_sedimnt.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/StratAer/aer_sedimnt.F90 (revision 2690) +++ LMDZ5/trunk/libf/phylmd/StratAer/aer_sedimnt.F90 (revision 2690) @@ -0,0 +1,130 @@ +SUBROUTINE AER_SEDIMNT(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) + +!**** *AER_SEDIMNT* - ROUTINE FOR PARAMETRIZATION OF AEROSOL SEDIMENTATION + +! Christoph Kleinschmitt +! based on the sedimentation scheme of +! Olivier Boucher & Jean-Jacques Morcrette +! (following the ice sedimentation scheme of Adrian Tompkins) + +!** INTERFACE. +! ---------- +! *AER_SEDIMNT* IS CALLED FROM *traccoag_mod*. + +!----------------------------------------------------------------------- + + USE phys_local_var_mod, ONLY: mdw, sfluxaer, DENSO4, f_r_wet, vsed_aer + USE dimphy, ONLY : klon,klev + USE infotrac + USE aerophys + USE YOMCST + +IMPLICIT NONE + +!----------------------------------------------------------------------- + + ! transfer variables when calling this routine + REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) + REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT):: tr_seri ! Concentration Traceur [U/KgA] + REAL,DIMENSION(klon,klev) :: dens_aer! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass + + ! local variables in sedimentation routine + INTEGER :: JL,JK,nb + REAL,DIMENSION(klon,klev) :: zvis ! dynamic viscosity of air [kg/(m*s)] + REAL,DIMENSION(klon,klev) :: zlair ! mean free path of air [m] + REAL :: ZRHO ! air density [kg/m^3] + REAL :: ZGDP ! =g/dp=1/(rho*dz) + REAL :: ZDTGDP ! =dt/(rho*dz) + REAL,DIMENSION(klon,nbtr_bin) :: ZSEDFLX ! sedimentation flux of tracer [U/(m^2*s)] + REAL,DIMENSION(nbtr_bin) :: ZAERONW ! tracer concentration at current time step [U/KgA] + REAL,DIMENSION(klon,nbtr_bin) :: ZAERONWM1! tracer concentration at preceding time step [U/KgA] + REAL,DIMENSION(klon,klev,nbtr_bin) :: ZVAER ! sedimentation velocity [m/s] + REAL,DIMENSION(nbtr_bin) :: ZSOLAERS ! sedimentation flux arriving from above [U/(m^2*s)] + REAL,DIMENSION(nbtr_bin) :: ZSOLAERB ! sedimentation flux leaving gridbox [U/(m^2*s)] + REAL,DIMENSION(klon,klev) :: m_sulf + +! dynamic viscosity of air (Pruppacher and Klett, 1978) [kg/(m*s)] +WHERE (t_seri.GE.273.15) + zvis=(1.718 + 0.0049*(t_seri-273.15))*1.E-5 + ELSEWHERE + zvis=(1.718 + 0.0049*(t_seri-273.15)-1.2E-05*(t_seri-273.15)**2)*1.E-5 +END WHERE + +! mean free path of air (Prupp. Klett) [m] +zlair(:,:) = 0.066 *(1.01325E+5/pplay(:,:))*(t_seri(:,:)/293.15)*1.E-06 + +!--initialisations of variables carried out from one layer to the next layer +!--actually not needed if (JK>1) test is on +DO JL=1,klon + DO nb=1,nbtr_bin + ZSEDFLX(JL,nb)=0.0 + ZAERONWM1(JL,nb)=0.0 + ENDDO +ENDDO + +!--from top to bottom (!) +DO JK=klev,1,-1 + DO JL=1,klon + DO nb=1,nbtr_bin + !--initialisations + ZSOLAERS(nb)=0.0 + ZSOLAERB(nb)=0.0 + ZGDP=RG/(paprs(JL,JK)-paprs(JL,JK+1)) + ZDTGDP=pdtphys*ZGDP + + ! source from above + IF (JK4) + REAL :: x ! molefraction of H2SO4 in the critical cluster + REAL Vbin(nbtr_bin) + REAL a_xm, b_xm, c_xm + REAL PDT, dt + REAL H2SO4_init + REAL ACTSO4(klon,klev) + REAL RRSI(nbtr_bin) + REAL nucl_rate + REAL cond_evap_rate + REAL evap_rate + REAL FL(nbtr_bin) + REAL ASO4(nbtr_bin) + REAL DNDR(nbtr_bin) + REAL H2SO4_sat(nbtr_bin) + + DO IK=1,nbtr_bin + Vbin(IK)=4.0*RPI*((mdw(IK)/2.)**3)/3.0 + ENDDO + + !coefficients for H2SO4 density parametrization used for nucleation if ntot<4 + a_xm = 0.7681724 + 1.*(2.1847140 + 1.*(7.1630022 + 1.*(-44.31447 + & + & 1.*(88.75606 + 1.*(-75.73729 + 1.*23.43228))))) + b_xm = 1.808225e-3 + 1.*(-9.294656e-3 + 1.*(-0.03742148 + 1.*(0.2565321 + & + & 1.*(-0.5362872 + 1.*(0.4857736 - 1.*0.1629592))))) + c_xm = -3.478524e-6 + 1.*(1.335867e-5 + 1.*(5.195706e-5 + 1.*(-3.717636e-4 + & + & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 ))))) + + ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap + CALL STRAACT(ACTSO4) + + ! compute particle radius in cm RRSI from diameter in m + DO it=1,nbtr_bin + RRSI(it)=mdw(it)/2.*100. + ENDDO + + DO ilon=1, klon + DO ilev=1, klev + ! only in the stratosphere + IF (is_strato(ilon,ilev)) THEN + ! initialize sulfur fluxes + sulf_nucl(ilon,ilev)=0.0 + sulf_cond_evap(ilon,ilev)=0.0 + H2SO4_init=tr_seri(ilon,ilev,id_H2SO4_strat) + ! adaptive timestep for nucleation and condensation + PDT=pdtphys + count_tstep=0 + DO while(PDT>0.0) + count_tstep=count_tstep+1 + if( count_tstep .GT. nbtstep ) exit + ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3) + rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) & + & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol + ! compute nucleation rate in kg(H2SO4)/kgA/s + CALL nucleation_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev),rh(ilon,ilev), & + & a_xm,b_xm,c_xm,nucl_rate,ntot,x) + ! compute cond/evap rate in kg(H2SO4)/kgA/s + CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & + & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & + & RRSI,Vbin,FL,ASO4,DNDR) + ! consider only condensation (positive FL) + DO IK=1,nbtr_bin + FL(IK)=MAX(FL(IK),0.) + ENDDO + ! compute total H2SO4 cond flux for all particles + cond_evap_rate=0.0 + DO IK=1, nbtr_bin + cond_evap_rate=cond_evap_rate+tr_seri(ilon,ilev,IK+nbtr_sulgas)*FL(IK)*mH2SO4mol + ENDDO + ! determine appropriate time step + dt=(H2SO4_init-H2SO4_sat(nbtr_bin))/float(nbtstep)/MAX(1.e-30, nucl_rate+cond_evap_rate) !cond_evap_rate pos. for cond. and neg. for evap. + IF (dt.LT.0.0) THEN + dt=PDT + ENDIF + dt=MIN(dt,PDT) + ! update H2SO4 concentration + tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-(nucl_rate+cond_evap_rate)*dt) + ! apply cond to bins + CALL cond_evap_part(dt,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:)) + ! apply nucl. to bins + CALL nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri(ilon,ilev,:)) + ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond) + sulf_cond_evap(ilon,ilev)=sulf_cond_evap(ilon,ilev)+mSatom/mH2SO4mol & + & *cond_evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*dt/pdtphys + sulf_nucl(ilon,ilev)=sulf_nucl(ilon,ilev)+mSatom/mH2SO4mol & + & *nucl_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*dt/pdtphys + ! update time step + PDT=PDT-dt + ENDDO + ! convert tr_seri(GASH2SO4) (in kg/kgA) to H2SO4 number density (in molecules/cm3) + rhoa=tr_seri(ilon,ilev,id_H2SO4_strat) & + & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol + ! compute cond/evap rate in kg(H2SO4)/kgA/s (now only evap for pdtphys) + CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & + & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & + & RRSI,Vbin,FL,ASO4,DNDR) + ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet + DO IK=1,nbtr_bin + FL(IK)=MAX(FL(IK)*pdtphys,0.-ASO4(IK))/pdtphys + ! consider only evap (negative FL) + FL(IK)=MIN(FL(IK),0.) + ENDDO + ! compute total H2SO4 evap flux for all particles + evap_rate=0.0 + DO IK=1, nbtr_bin + evap_rate=evap_rate+tr_seri(ilon,ilev,IK+nbtr_sulgas)*FL(IK)*mH2SO4mol + ENDDO + ! update H2SO4 concentration after evap + tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-evap_rate*pdtphys) + ! apply evap to bins + CALL cond_evap_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:)) + ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond) + sulf_cond_evap(ilon,ilev)=sulf_cond_evap(ilon,ilev)+mSatom/mH2SO4mol & + & *evap_rate*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG + ENDIF + ENDDO + ENDDO + + IF (MINVAL(tr_seri(:,:,:)).LT.0.0) THEN + DO ilon=1, klon + DO ilev=1, klev + DO IK=1, nbtr + IF (tr_seri(ilon,ilev,IK).LT.0.0) THEN + PRINT *, 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,IK), ilon, ilev, IK + ENDIF + ENDDO + ENDDO + ENDDO + ENDIF + +END SUBROUTINE micphy_tstep Index: LMDZ5/trunk/libf/phylmd/StratAer/miecalc_aer.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/StratAer/miecalc_aer.F90 (revision 2690) +++ LMDZ5/trunk/libf/phylmd/StratAer/miecalc_aer.F90 (revision 2690) @@ -0,0 +1,564 @@ +SUBROUTINE MIECALC_AER(tau_strat, piz_strat, cg_strat, tau_strat_wave, tau_lw_abs_rrtm, paprs, debut) + +!-------Mie computations for a size distribution +! of homogeneous spheres. +! +!========================================================== +!--Ref : Toon and Ackerman, Applied Optics, 1981 +! Stephens, CSIRO, 1979 +! Attention : surdimensionement des tableaux +! to be compiled with double precision option (-r8 on Sun) +! AUTHOR: Olivier Boucher, Christoph Kleinschmitt +!-------SIZE distribution properties---------------- +!--sigma_g : geometric standard deviation +!--r_0 : geometric number mean radius (um)/modal radius +!--Ntot : total concentration in m-3 + + USE phys_local_var_mod, ONLY: tr_seri, mdw, alpha_bin, piz_bin, cg_bin + USE aerophys + USE aero_mod + USE infotrac, ONLY : nbtr, nbtr_bin, nbtr_sulgas, id_SO2_strat + USE dimphy + USE YOMCST , ONLY : RG, RPI + USE mod_phys_lmdz_para, only: gather, scatter, bcast + USE mod_grid_phy_lmdz, ONLY : klon_glo + USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root + + IMPLICIT NONE + +! Variable input + LOGICAL,INTENT(IN) :: debut ! le flag de l'initialisation de la physique + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) + +! Stratospheric aerosols optical properties + REAL, DIMENSION(klon,klev,nbands_sw_rrtm) :: tau_strat, piz_strat, cg_strat + REAL, DIMENSION(klon,klev,nwave_sw+nwave_lw) :: tau_strat_wave + REAL, DIMENSION(klon,klev,nbands_lw_rrtm) :: tau_lw_abs_rrtm + +!! REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA] + +! local variables + REAL Ntot + PARAMETER (Ntot=1.0) + LOGICAL, PARAMETER :: refr_ind_interpol = .TRUE. ! set interpolation of refractive index + REAL r_0 ! aerosol particle radius [m] + INTEGER bin_number, ilon, ilev + REAL masse,volume,surface + REAL rmin, rmax !----integral bounds in m + +!------------------------------------- + + COMPLEX m !----refractive index m=n_r-i*n_i + INTEGER Nmax,Nstart !--number of iterations + COMPLEX k2, k3, z1, z2 + COMPLEX u1,u5,u6,u8 + COMPLEX a(1:21000), b(1:21000) + COMPLEX I + INTEGER n !--loop index + REAL nnn + COMPLEX nn + REAL Q_ext, Q_abs, Q_sca, g, omega !--parameters for radius r + REAL x !--size parameter + REAL r, r_lower, r_upper !--radius + REAL sigma_sca, sigma_ext, sigma_abs + REAL omegatot, gtot !--averaged parameters + COMPLEX ksiz2(-1:21000), psiz2(1:21000) + COMPLEX nu1z1(1:21010), nu1z2(1:21010) + COMPLEX nu3z2(0:21000) + REAL number, deltar + INTEGER bin, Nbin, it + PARAMETER (Nbin=10) + +!--has to be consistent with dataset + INTEGER nb_lambda_h2so4 + PARAMETER (nb_lambda_h2so4=62) !61, 235 + +!---wavelengths STREAMER + INTEGER Nwv, NwvmaxSW + PARAMETER (NwvmaxSW=24) + REAL lambda(1:NwvmaxSW+1) + DATA lambda/0.28E-6, 0.30E-6, 0.33E-6, 0.36E-6, 0.40E-6, & + 0.44E-6, 0.48E-6, 0.52E-6, 0.57E-6, 0.64E-6, & + 0.69E-6, 0.75E-6, 0.78E-6, 0.87E-6, 1.00E-6, & + 1.10E-6, 1.19E-6, 1.28E-6, 1.53E-6, 1.64E-6, & + 2.13E-6, 2.38E-6, 2.91E-6, 3.42E-6, 4.00E-6/ + +!---wavelengths de references +!---be careful here the 5th wavelength is 1020 nm + INTEGER nb + REAL lambda_ref(nwave_sw+nwave_lw) + DATA lambda_ref /0.443E-6,0.550E-6,0.670E-6, & + 0.765E-6,1.020E-6,10.E-6/ + +!--LW + INTEGER NwvmaxLW + PARAMETER (NwvmaxLW=500) + REAL Tb, hh, cc, kb + PARAMETER (Tb=220.0, hh=6.62607e-34) + PARAMETER (cc=2.99792e8, kb=1.38065e-23) + +!---TOA fluxes - Streamer Cs + REAL weight(1:NwvmaxSW), weightLW +!c DATA weight/0.839920E1, 0.231208E2, 0.322393E2, 0.465058E2, +!c . 0.678199E2, 0.798964E2, 0.771359E2, 0.888472E2, +!c . 0.115281E3, 0.727565E2, 0.816992E2, 0.336172E2, +!c . 0.914603E2, 0.112706E3, 0.658840E2, 0.524470E2, +!c . 0.391067E2, 0.883864E2, 0.276672E2, 0.681812E2, +!c . 0.190966E2, 0.250766E2, 0.128704E2, 0.698720E1/ +!---TOA fluxes - Tad + DATA weight/ 4.20, 11.56, 16.12, 23.25, 33.91, 39.95, 38.57, & + 44.42, 57.64, 29.36, 47.87, 16.81, 45.74, 56.35, & + 32.94, 26.22, 19.55, 44.19, 13.83, 34.09, 9.55, & + 12.54, 6.44, 3.49/ +!C---BOA fluxes - Tad +!c DATA weight/ 0.01, 4.05, 9.51, 15.99, 26.07, 33.10, 33.07, +!c . 39.91, 52.67, 27.89, 43.60, 13.67, 42.22, 40.12, +!c . 32.70, 14.44, 19.48, 14.23, 13.43, 16.42, 8.33, +!c . 0.95, 0.65, 2.76/ + + REAL lambda_int(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW), ll + REAL dlambda_int(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW), dl + + REAL n_r(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW) + REAL n_i(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW) + + REAL ilambda, ilambda_prev, ilambda_max, ilambda_min + REAL n_r_h2so4, n_i_h2so4 + REAL n_r_h2so4_prev, n_i_h2so4_prev + + REAL final_a(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW) + REAL final_g(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW) + REAL final_w(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW) + + INTEGER band, bandSW, bandLW + +!---wavelengths SW RRTM + REAL wv_rrtm_SW(nbands_sw_rrtm+1) + DATA wv_rrtm_SW/ 0.185E-6, 0.25E-6, 0.44E-6, 0.69E-6, & + 1.19E-6, 2.38E-6, 4.00E-6/ + +!---wavenumbers and wavelengths LW RRTM + REAL wn_rrtm(nbands_lw_rrtm+1), wv_rrtm(nbands_lw_rrtm+1) + DATA wn_rrtm/ 10., 250., 500., 630., 700., 820., & + 980., 1080., 1180., 1390., 1480., 1800., & + 2080., 2250., 2380., 2600., 3000./ + +!--GCM results + REAL gcm_a(nbands_sw_rrtm+nbands_lw_rrtm) + REAL gcm_g(nbands_sw_rrtm+nbands_lw_rrtm) + REAL gcm_w(nbands_sw_rrtm+nbands_lw_rrtm) + REAL gcm_weight_a(nbands_sw_rrtm+nbands_lw_rrtm) + REAL gcm_weight_g(nbands_sw_rrtm+nbands_lw_rrtm) + REAL gcm_weight_w(nbands_sw_rrtm+nbands_lw_rrtm) + + REAL ss_a(nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw) + REAL ss_w(nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw) + REAL ss_g(nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw) + + REAL wavenumber + REAL, DIMENSION(nb_lambda_h2so4,4) :: ref_ind + +!--------------------------------------------------------- + + IF (debut) THEN + !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994) + mdw(1)=mdwmin + IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio + mdw(2)=mdw(1)*2.**(1./3.) + DO it=3, nbtr_bin + mdw(it)=mdw(it-1)*V_rat**(1./3.) + ENDDO + ELSE + DO it=2, nbtr_bin + mdw(it)=mdw(it-1)*V_rat**(1./3.) + ENDDO + ENDIF + PRINT *,'init mdw=', mdw + +!$OMP MASTER + +! CALL gather(tr_seri, tr_seri_glo) +! IF (is_mpi_root) THEN +! IF (MAXVAL(tr_seri_glo).LT.1e-30) THEN + + !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K + DO bin_number=1, nbtr_bin + r_0=(dens_aer_dry/dens_aer_ref/0.75)**(1./3.)*mdw(bin_number)/2. + !--integral boundaries set to bin boundaries + rmin=r_0/sqrt(V_rat**(1./3.)) + rmax=r_0*sqrt(V_rat**(1./3.)) + + !--set up SW + DO Nwv=1, NwvmaxSW + lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2. + ENDDO + + DO nb=1, nwave_sw+nwave_lw + lambda_int(NwvmaxSW+nb)=lambda_ref(nb) + ENDDO + + !--set up LW + !--conversion wavenumber in cm-1 to wavelength in m + DO Nwv=1, nbands_lw_rrtm+1 + wv_rrtm(Nwv)=10000./wn_rrtm(Nwv)*1.e-6 + ENDDO + + DO Nwv=1, NwvmaxLW + lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= & + exp( log(wv_rrtm(1))+float(Nwv-1)/float(NwvmaxLW-1)* & + (log(wv_rrtm(nbands_lw_rrtm+1))-log(wv_rrtm(1))) ) + ENDDO + +!--computing the dlambdas + Nwv=1 + dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= & + & lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)- & + & lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv+1) + DO Nwv=2, NwvmaxLW-1 + dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= & + & (lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv-1)- & + & lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv+1))/2. + ENDDO + Nwv=NwvmaxLW + dlambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)= & + & lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv-1)- & + & lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv) + + OPEN (unit=11,file='h2so4_0.75_300.00_hummel_1988_p_q.dat') + + IF (refr_ind_interpol) THEN + DO nb=1,nb_lambda_h2so4 + READ(11,*) ref_ind(nb,:) + ENDDO + ilambda_max=ref_ind(1,2)/1.e6 !--in m + ilambda_min=ref_ind(nb_lambda_h2so4,2)/1.e6 !--in m + DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW + IF (lambda_int(Nwv).GT.ilambda_max) THEN + !for lambda out of data range, take boundary values + n_r(Nwv)=ref_ind(1,3) + n_i(Nwv)=ref_ind(1,4) + ELSEIF (lambda_int(Nwv).LE.ilambda_min) THEN + n_r(Nwv)=ref_ind(nb_lambda_h2so4,3) + n_i(Nwv)=ref_ind(nb_lambda_h2so4,4) + ELSE + DO nb=2,nb_lambda_h2so4 + ilambda=ref_ind(nb,2)/1.e6 + ilambda_prev=ref_ind(nb-1,2)/1.e6 + n_r_h2so4=ref_ind(nb,3) + n_r_h2so4_prev=ref_ind(nb-1,3) + n_i_h2so4=ref_ind(nb,4) + n_i_h2so4_prev=ref_ind(nb-1,4) + IF (lambda_int(Nwv).GT.ilambda.AND. & + lambda_int(Nwv).LE.ilambda_prev) THEN !--- linear interpolation + n_r(Nwv)=n_r_h2so4+(lambda_int(Nwv)-ilambda)/ & + (ilambda_prev-ilambda)*(n_r_h2so4_prev-n_r_h2so4) + n_i(Nwv)=n_i_h2so4+(lambda_int(Nwv)-ilambda)/ & + (ilambda_prev-ilambda)*(n_i_h2so4_prev-n_i_h2so4) + ENDIF + ENDDO + ENDIF + ENDDO + ELSE + DO nb=1,nb_lambda_h2so4 + READ(11,*) wavenumber, ilambda, n_r_h2so4, n_i_h2so4 + ilambda=1000.*ilambda !wavelength in nm + DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW + IF (ilambda/1.e9.GT.lambda_int(Nwv)) THEN !--- step function + n_r(Nwv)=n_r_h2so4 + n_i(Nwv)=n_i_h2so4 + ENDIF + ENDDO + ENDDO + ENDIF + CLOSE(11) + + !---Loop on wavelengths + DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW + + m=CMPLX(n_r(Nwv),-n_i(Nwv)) + + I=CMPLX(0.,1.) + + sigma_sca=0.0 + sigma_ext=0.0 + sigma_abs=0.0 + gtot=0.0 + omegatot=0.0 + masse = 0.0 + volume=0.0 + surface=0.0 + + DO bin=1, Nbin !---loop on size bins + + r_lower=exp(log(rmin)+FLOAT(bin-1)/FLOAT(Nbin)*(log(rmax)-log(rmin))) + r_upper=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin))) + r=sqrt(r_lower*r_upper) + x=2.*RPI*r/lambda_int(Nwv) + deltar=r_upper-r_lower + + number=Ntot*deltar/(rmax-rmin) !dN/dr constant over tracer bin +! masse=masse +4./3.*RPI*(r**3)*number*deltar*ropx*1.E3 !--g/m3 + volume=volume+4./3.*RPI*(r**3)*number*deltar + surface=surface+4.*RPI*r**2*number*deltar + + k2=m + k3=CMPLX(1.0,0.0) + + z2=CMPLX(x,0.0) + z1=m*z2 + + IF (0.0.LE.x.AND.x.LE.8.) THEN + Nmax=INT(x+4*x**(1./3.)+1.)+2 + ELSEIF (8..LT.x.AND.x.LT.4200.) THEN + Nmax=INT(x+4.05*x**(1./3.)+2.)+1 + ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN + Nmax=INT(x+4*x**(1./3.)+2.)+1 + ELSE + PRINT *, 'x out of bound, x=', x + STOP + ENDIF + + Nstart=Nmax+10 + + !-----------loop for nu1z1, nu1z2 + + nu1z1(Nstart)=CMPLX(0.0,0.0) + nu1z2(Nstart)=CMPLX(0.0,0.0) + DO n=Nstart-1, 1 , -1 + nn=CMPLX(FLOAT(n),0.0) + nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) ) + nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) ) + ENDDO + + !------------loop for nu3z2 + + nu3z2(0)=-I + DO n=1, Nmax + nn=CMPLX(FLOAT(n),0.0) + nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) ) + ENDDO + + !-----------loop for psiz2 and ksiz2 (z2) + ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2)) + ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2)) + DO n=1,Nmax + nn=CMPLX(FLOAT(n),0.0) + ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2) + psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0) + ENDDO + + !-----------loop for a(n) and b(n) + + DO n=1, Nmax + u1=k3*nu1z1(n) - k2*nu1z2(n) + u5=k3*nu1z1(n) - k2*nu3z2(n) + u6=k2*nu1z1(n) - k3*nu1z2(n) + u8=k2*nu1z1(n) - k3*nu3z2(n) + a(n)=psiz2(n)/ksiz2(n) * u1/u5 + b(n)=psiz2(n)/ksiz2(n) * u6/u8 + ENDDO + + !-----------------final loop-------------- + Q_ext=0.0 + Q_sca=0.0 + g=0.0 + DO n=Nmax-1,1,-1 + nnn=FLOAT(n) + Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) ) + Q_sca=Q_sca+ (2.*nnn+1.) * & + REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) ) + g=g + nnn*(nnn+2.)/(nnn+1.) * & + REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) ) + & + (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n))) + ENDDO + Q_ext=2./x**2 * Q_ext + Q_sca=2./x**2 * Q_sca + Q_abs=Q_ext-Q_sca + IF (AIMAG(m).EQ.0.0) Q_abs=0.0 + omega=Q_sca/Q_ext + g=g*4./x**2/Q_sca + + sigma_sca=sigma_sca+r**2*Q_sca*number + sigma_abs=sigma_abs+r**2*Q_abs*number + sigma_ext=sigma_ext+r**2*Q_ext*number + omegatot=omegatot+r**2*Q_ext*omega*number + gtot =gtot+r**2*Q_sca*g*number + + ENDDO !---bin + !------------------------------------------------------------------ + + sigma_sca=RPI*sigma_sca + sigma_abs=RPI*sigma_abs + sigma_ext=RPI*sigma_ext + gtot=RPI*gtot/sigma_sca + omegatot=RPI*omegatot/sigma_ext + + final_g(Nwv)=gtot + final_w(Nwv)=omegatot +! final_a(Nwv)=sigma_ext/masse + final_a(Nwv)=sigma_ext !extinction/absorption cross section per particle + + ENDDO !--loop on wavelength + + + !---averaging over LMDZ wavebands + + DO band=1, nbands_sw_rrtm+nbands_lw_rrtm + gcm_a(band)=0.0 + gcm_g(band)=0.0 + gcm_w(band)=0.0 + gcm_weight_a(band)=0.0 + gcm_weight_g(band)=0.0 + gcm_weight_w(band)=0.0 + ENDDO + + !---band 1 is now in the UV, so we take the first wavelength as being representative + DO Nwv=1,1 + bandSW=1 + gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv) + gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv) + gcm_w(bandSW)=gcm_w(bandSW)+final_w(Nwv)*final_a(Nwv)*weight(Nwv) + gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+final_a(Nwv)*weight(Nwv) + gcm_g(bandSW)=gcm_g(bandSW)+final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv) + gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+final_a(Nwv)*final_w(Nwv)*weight(Nwv) + ENDDO + + DO Nwv=1,NwvmaxSW + + IF (lambda_int(Nwv).LE.wv_rrtm_SW(3)) THEN !--RRTM spectral interval 2 + bandSW=2 + ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(4)) THEN !--RRTM spectral interval 3 + bandSW=3 + ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(5)) THEN !--RRTM spectral interval 4 + bandSW=4 + ELSEIF (lambda_int(Nwv).LE.wv_rrtm_SW(6)) THEN !--RRTM spectral interval 5 + bandSW=5 + ELSE !--RRTM spectral interval 6 + bandSW=6 + ENDIF + + gcm_a(bandSW)=gcm_a(bandSW)+final_a(Nwv)*weight(Nwv) + gcm_weight_a(bandSW)=gcm_weight_a(bandSW)+weight(Nwv) + gcm_w(bandSW)=gcm_w(bandSW)+final_w(Nwv)*final_a(Nwv)*weight(Nwv) + gcm_weight_w(bandSW)=gcm_weight_w(bandSW)+final_a(Nwv)*weight(Nwv) + gcm_g(bandSW)=gcm_g(bandSW)+final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv) + gcm_weight_g(bandSW)=gcm_weight_g(bandSW)+final_a(Nwv)*final_w(Nwv)*weight(Nwv) + + ENDDO + + DO Nwv=NwvmaxSW+nwave_sw+nwave_lw+1,NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW + ll=lambda_int(Nwv) + dl=dlambda_int(Nwv) + weightLW=1./ll**5./(exp(hh*cc/kb/Tb/ll)-1.)*dl + bandLW=1 !--default value starting from the highest lambda + DO band=2, nbands_lw_rrtm + IF (ll.LT.wv_rrtm(band)) THEN !--as long as + bandLW=band + ENDIF + ENDDO + gcm_a(nbands_sw_rrtm+bandLW)=gcm_a(nbands_sw_rrtm+bandLW)+final_a(Nwv)* & + (1.-final_w(Nwv))*weightLW + gcm_weight_a(nbands_sw_rrtm+bandLW)=gcm_weight_a(nbands_sw_rrtm+bandLW)+weightLW + + gcm_w(nbands_sw_rrtm+bandLW)=gcm_w(nbands_sw_rrtm+bandLW)+final_w(Nwv)* & + final_a(Nwv)*weightLW + gcm_weight_w(nbands_sw_rrtm+bandLW)=gcm_weight_w(nbands_sw_rrtm+bandLW)+ & + final_a(Nwv)*weightLW + + gcm_g(nbands_sw_rrtm+bandLW)=gcm_g(nbands_sw_rrtm+bandLW)+final_g(Nwv)* & + final_a(Nwv)*final_w(Nwv)*weightLW + gcm_weight_g(nbands_sw_rrtm+bandLW)=gcm_weight_g(nbands_sw_rrtm+bandLW)+ & + final_a(Nwv)*final_w(Nwv)*weightLW + ENDDO + + DO band=1, nbands_sw_rrtm+nbands_lw_rrtm + gcm_a(band)=gcm_a(band)/gcm_weight_a(band) + gcm_w(band)=gcm_w(band)/gcm_weight_w(band) + gcm_g(band)=gcm_g(band)/gcm_weight_g(band) + ss_a(band)=gcm_a(band) + ss_w(band)=gcm_w(band) + ss_g(band)=gcm_g(band) + ENDDO + + DO nb=1, nwave_sw+nwave_lw + ss_a(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_a(NwvmaxSW+nb) + ss_w(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_w(NwvmaxSW+nb) + ss_g(nbands_sw_rrtm+nbands_lw_rrtm+nb)=final_g(NwvmaxSW+nb) + ENDDO + + DO nb=1,nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw + alpha_bin(nb,bin_number)=ss_a(nb) !extinction/absorption cross section per particle + piz_bin(nb,bin_number)=ss_w(nb) + cg_bin(nb,bin_number)=ss_g(nb) + + ENDDO + + ENDDO !loop over tracer bins + +! ENDIF !mpi_root + +!$OMP END MASTER + CALL bcast(alpha_bin) + CALL bcast(piz_bin) + CALL bcast(cg_bin) +!$OMP BARRIER + +! CALL scatter(alpha_bin, alpha_bin) +! CALL scatter(piz_bin, piz_bin) +! CALL scatter(cg_bin, cg_bin) + + !set to default values at first time step (tr_seri still zero) + tau_strat(:,:,:)=1.e-15 + piz_strat(:,:,:)=1.0 + cg_strat(:,:,:)=0.0 + tau_lw_abs_rrtm(:,:,:)=1.e-15 + tau_strat_wave(:,:,:)=1.e-15 + + ELSE + +! CALL gather(tr_seri, tr_seri_glo) + + !--compute optical properties of actual size distribution (from tr_seri) + DO ilon=1,klon + DO ilev=1, klev + DO nb=1,nbands_sw_rrtm + tau_strat(ilon,ilev,nb)=0.0 + DO bin_number=1, nbtr_bin + tau_strat(ilon,ilev,nb)=tau_strat(ilon,ilev,nb)+alpha_bin(nb,bin_number) & + *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG + ENDDO + + piz_strat(ilon,ilev,nb)=0.0 + DO bin_number=1, nbtr_bin + piz_strat(ilon,ilev,nb)=piz_strat(ilon,ilev,nb)+piz_bin(nb,bin_number)*alpha_bin(nb,bin_number) & + *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG + ENDDO + piz_strat(ilon,ilev,nb)=piz_strat(ilon,ilev,nb)/MAX(tau_strat(ilon,ilev,nb),1.e-15) + + cg_strat(ilon,ilev,nb)=0.0 + DO bin_number=1, nbtr_bin + cg_strat(ilon,ilev,nb)=cg_strat(ilon,ilev,nb)+cg_bin(nb,bin_number)*piz_bin(nb,bin_number)*alpha_bin(nb,bin_number) & + *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG + ENDDO + cg_strat(ilon,ilev,nb)=cg_strat(ilon,ilev,nb)/MAX(tau_strat(ilon,ilev,nb)*piz_strat(ilon,ilev,nb),1.e-15) + ENDDO + DO nb=1,nbands_lw_rrtm + tau_lw_abs_rrtm(ilon,ilev,nb)=0.0 + DO bin_number=1, nbtr_bin + tau_lw_abs_rrtm(ilon,ilev,nb)=tau_lw_abs_rrtm(ilon,ilev,nb)+alpha_bin(nbands_sw_rrtm+nb,bin_number) & + *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG + ENDDO + ENDDO + DO nb=1,nwave_sw+nwave_lw + tau_strat_wave(ilon,ilev,nb)=0.0 + DO bin_number=1, nbtr_bin + tau_strat_wave(ilon,ilev,nb)=tau_strat_wave(ilon,ilev,nb)+alpha_bin(nbands_sw_rrtm+nbands_lw_rrtm+nb,bin_number) & + *tr_seri(ilon,ilev,bin_number+nbtr_sulgas)*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG + ENDDO + ENDDO + ENDDO + ENDDO + + ENDIF !debut + +END SUBROUTINE MIECALC_AER Index: LMDZ5/trunk/libf/phylmd/StratAer/minmaxsimple.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/StratAer/minmaxsimple.F90 (revision 2690) +++ LMDZ5/trunk/libf/phylmd/StratAer/minmaxsimple.F90 (revision 2690) @@ -0,0 +1,23 @@ +! +! $Id: minmaxsimple.F90 1910 2013-11-29 08:40:25Z fairhead $ +! +SUBROUTINE minmaxsimple(zq,qmin,qmax,comment) + USE dimphy + IMPLICIT NONE + +! Entrees + REAL,DIMENSION(klon,klev), INTENT(IN) :: zq + REAL,INTENT(IN) :: qmin,qmax + CHARACTER(LEN=*),INTENT(IN) :: comment + +! Local + REAL zmin, zmax + + zmin=MINVAL(zq) + zmax=MAXVAL(zq) + PRINT *, "qmin qmax=", zmin, zmax, comment +! IF (zmin.LT.qmin.OR.zmax.GT.qmax) THEN +! WRITE(*,*) "qmin qmax=", zmin, zmax, comment +! ENDIF + +END SUBROUTINE minmaxsimple Index: LMDZ5/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90 (revision 2690) +++ LMDZ5/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90 (revision 2690) @@ -0,0 +1,278 @@ +MODULE nucleation_tstep_mod + +CONTAINS + +SUBROUTINE nucleation_rate(rhoa,t_seri,pplay,rh,a_xm,b_xm,c_xm,nucl_rate,ntot,x) + + USE aerophys + USE infotrac + USE YOMCST + + IMPLICIT NONE + + ! input variables + REAL rhoa !H2SO4 number density [molecules/cm3] + REAL t_seri + REAL pplay + REAL rh + REAL a_xm, b_xm, c_xm + + ! output variables + REAL nucl_rate + REAL ntot ! total number of molecules in the critical cluster + REAL x ! molefraction of H2SO4 in the critical cluster + + ! local variables + REAL, PARAMETER :: k_B=1.3806E-23 ! Boltzmann constant [J/K] + REAL :: jnuc !nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s) + REAL :: rc !radius of the critical cluster in nm + REAL VH2SO4mol + + ! call nucleation routine + CALL binapara(t_seri,rh,rhoa,jnuc,x,ntot,rc) + + IF (ntot < 4.0) THEN + !set jnuc to collision rate of two H2SO4 molecules (following personal communication of Ulrike Niemeier and Hanna Vehkamäki) + VH2SO4mol=mH2SO4mol/(1.e-3*(a_xm+t_seri*(b_xm+t_seri*c_xm))) !cm3 + jnuc = rhoa**2. *(3./4.*RPI)**(1./6.) *(12.*k_B*t_seri/mH2SO4mol)**0.5 & + & *100.*(2.*VH2SO4mol**(1./3.))**2. !1/(cm3s) + ntot=2.0 + x=1.0 + ENDIF + + ! convert jnuc from particles/cm3/s to kg(H2SO4)/kgA/s + nucl_rate=jnuc*ntot*x*mH2SO4mol/(pplay/t_seri/RD/1.E6) + +END SUBROUTINE nucleation_rate + +!-------------------------------------------------------------------------------------------------- + +SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri) + + USE aerophys + USE infotrac + USE YOMCST + + IMPLICIT NONE + + ! input variables + REAL nucl_rate + REAL ntot ! total number of molecules in the critical cluster + REAL x ! molefraction of H2SO4 in the critical cluster + REAL dt + REAL Vbin(nbtr_bin) + + ! output variables + REAL tr_seri(nbtr) + + ! local variables + INTEGER k + REAL Vnew + REAL ff(nbtr_bin) + + ! dry volume of nucleated particle (at reference temperature) + Vnew=mH2SO4mol*ntot*x/dens_aer_dry + + ! compute distribution factor for particles of intermediate size (from Jacobson 1994, equation 13) + ff(:)=0.0 + DO k=1, nbtr_bin + ! CK 20160531: bug fix for first bin + IF (k.LE.(nbtr_bin-1).AND.Vbin(k).LE.Vnew.AND.Vnew.LT.Vbin(k+1)) THEN + ff(k)= Vbin(k)/Vnew*(Vbin(k+1)-Vnew)/(Vbin(k+1)-Vbin(k)) + ENDIF + IF (k.EQ.1.AND.Vnew.LE.Vbin(k)) THEN + ff(k)= 1. + ENDIF + IF (k.GT.1.AND.Vbin(k-1).LT.Vnew.AND.Vnew.LT.Vbin(k)) THEN + ff(k)= 1.-ff(k-1) + ENDIF + IF (k.EQ.nbtr_bin.AND.Vnew.GE.Vbin(k)) THEN + ff(k)= 1. + ENDIF + ENDDO + ! correction of ff for volume conservation + DO k=1, nbtr_bin + ff(k)=ff(k)*Vnew/Vbin(k) + ENDDO + + ! add nucleated H2SO4 to corresponding aerosol bin + DO k=1, nbtr_bin + tr_seri(k+nbtr_sulgas)=tr_seri(k+nbtr_sulgas)+nucl_rate*dt/(mH2SO4mol*ntot*x)*ff(k) + ENDDO + +END SUBROUTINE nucleation_part + +!--------------------------------------------------------------------------------------------------- + +SUBROUTINE binapara(pt,prh,rhoa_in,jnuc,x,ntot,rc) + + + ! Fortran 90 subroutine binapara + ! + ! Calculates parametrized values of nucleation rate, + ! mole fraction of sulphuric acid + ! total number of particles, and the radius of the critical cluster + ! in H2O-H2SO4 system IF temperature, saturatio ratio of water and + ! sulfuric acid concentration are given. + ! + ! Copyright (C) 2002 Hanna Vehkamäki + ! + ! Division of Atmospheric Sciences + ! Department of Physical Sciences + ! P.O. Box 64 + ! FIN-00014 University of Helsinki + ! Finland + ! + ! hanna.vehkamaki@helsinki.fi + + IMPLICIT NONE + + REAL :: pt,t !temperature in K (190.15-300.15K) + REAL :: prh,rh !saturatio ratio of water (0.0001-1) + REAL,intent(in) :: rhoa_in !sulfuric acid concentration in 1/cm3 (10^4-10^11 1/cm3) + REAL,intent(out) :: jnuc !nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s) + REAL,intent(out) :: ntot !total number of molecules in the critical cluster (ntot>4) + REAL,intent(out) :: x ! molefraction of H2SO4 in the critical cluster + REAL,intent(out) :: rc !radius of the critical cluster in nm + REAL :: rhotres ! treshold concentration of h2so4 (1/cm^3) + ! which produces nucleation rate 1/(cm^3 s) as a function of rh and t + REAL rhoa + +! CK: use intermediate variables to avoid overwriting + t=pt + rh=prh + rhoa=rhoa_in + + IF (t < 190.15) THEN +! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature < 190.15 K, using 190.15 K (T=',t,'K)' + t=190.15 + ENDIF + + IF (t > 300.15) THEN +! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature > 300.15 K, using 300.15 K' + t=300.15 + ENDIF + + IF (rh < 0.0001) THEN +! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water < 0.0001, using 0.0001' + rh=0.0001 + ENDIF + + IF (rh > 1.0) THEN +! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water > 1 using 1' +! print *, 'rh=',rh + rh=1.0 + ENDIF + + IF (rhoa < 1.e4) THEN +! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3' + rhoa=1.e4 + ENDIF + + IF (rhoa > 1.e11) THEN +! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration > 1e11 1/cm3, using 1e11 1/cm3' + rhoa=1.e11 + ENDIF + + x= 0.7409967177282139 - 0.002663785665140117*t + 0.002010478847383187*Log(rh) & + & - 0.0001832894131464668*t*Log(rh) + 0.001574072538464286*Log(rh)**2 & + & - 0.00001790589121766952*t*Log(rh)**2 + 0.0001844027436573778*Log(rh)**3 & + & - 1.503452308794887e-6*t*Log(rh)**3 - 0.003499978417957668*Log(rhoa) & + & + 0.0000504021689382576*t*Log(rhoa) + + jnuc= 0.1430901615568665 + 2.219563673425199*t - 0.02739106114964264*t**2 + & + & 0.00007228107239317088*t**3 + 5.91822263375044/x + & + & 0.1174886643003278*Log(rh) + 0.4625315047693772*t*Log(rh) - & + & 0.01180591129059253*t**2*Log(rh) + & + & 0.0000404196487152575*t**3*Log(rh) + (15.79628615047088*Log(rh))/x - & + & 0.215553951893509*Log(rh)**2 - 0.0810269192332194*t*Log(rh)**2 + & + & 0.001435808434184642*t**2*Log(rh)**2 - & + & 4.775796947178588e-6*t**3*Log(rh)**2 - & + & (2.912974063702185*Log(rh)**2)/x - 3.588557942822751*Log(rh)**3 + & + & 0.04950795302831703*t*Log(rh)**3 - & + & 0.0002138195118737068*t**2*Log(rh)**3 + & + & 3.108005107949533e-7*t**3*Log(rh)**3 - & + & (0.02933332747098296*Log(rh)**3)/x + & + & 1.145983818561277*Log(rhoa) - & + & 0.6007956227856778*t*Log(rhoa) + & + & 0.00864244733283759*t**2*Log(rhoa) - & + & 0.00002289467254710888*t**3*Log(rhoa) - & + & (8.44984513869014*Log(rhoa))/x + & + & 2.158548369286559*Log(rh)*Log(rhoa) + & + & 0.0808121412840917*t*Log(rh)*Log(rhoa) - & + & 0.0004073815255395214*t**2*Log(rh)*Log(rhoa) - & + & 4.019572560156515e-7*t**3*Log(rh)*Log(rhoa) + & + & (0.7213255852557236*Log(rh)*Log(rhoa))/x + & + & 1.62409850488771*Log(rh)**2*Log(rhoa) - & + & 0.01601062035325362*t*Log(rh)**2*Log(rhoa) + & + & 0.00003771238979714162*t**2*Log(rh)**2*Log(rhoa) + & + & 3.217942606371182e-8*t**3*Log(rh)**2*Log(rhoa) - & + & (0.01132550810022116*Log(rh)**2*Log(rhoa))/x + & + & 9.71681713056504*Log(rhoa)**2 - & + & 0.1150478558347306*t*Log(rhoa)**2 + & + & 0.0001570982486038294*t**2*Log(rhoa)**2 + & + & 4.009144680125015e-7*t**3*Log(rhoa)**2 + & + & (0.7118597859976135*Log(rhoa)**2)/x - & + & 1.056105824379897*Log(rh)*Log(rhoa)**2 + & + & 0.00903377584628419*t*Log(rh)*Log(rhoa)**2 - & + & 0.00001984167387090606*t**2*Log(rh)*Log(rhoa)**2 + & + & 2.460478196482179e-8*t**3*Log(rh)*Log(rhoa)**2 - & + & (0.05790872906645181*Log(rh)*Log(rhoa)**2)/x - & + & 0.1487119673397459*Log(rhoa)**3 + & + & 0.002835082097822667*t*Log(rhoa)**3 - & + & 9.24618825471694e-6*t**2*Log(rhoa)**3 + & + & 5.004267665960894e-9*t**3*Log(rhoa)**3 - & + & (0.01270805101481648*Log(rhoa)**3)/x + jnuc=exp(jnuc) !1/(cm3s) + + ntot =-0.002954125078716302 - 0.0976834264241286*t + 0.001024847927067835*t**2 - 2.186459697726116e-6*t**3 - & + & 0.1017165718716887/x - 0.002050640345231486*Log(rh) - 0.007585041382707174*t*Log(rh) + & + & 0.0001926539658089536*t**2*Log(rh) - 6.70429719683894e-7*t**3*Log(rh) - & + & (0.2557744774673163*Log(rh))/x + 0.003223076552477191*Log(rh)**2 + 0.000852636632240633*t*Log(rh)**2 - & + & 0.00001547571354871789*t**2*Log(rh)**2 + 5.666608424980593e-8*t**3*Log(rh)**2 + & + & (0.03384437400744206*Log(rh)**2)/x + 0.04743226764572505*Log(rh)**3 - & + & 0.0006251042204583412*t*Log(rh)**3 + 2.650663328519478e-6*t**2*Log(rh)**3 - & + & 3.674710848763778e-9*t**3*Log(rh)**3 - (0.0002672510825259393*Log(rh)**3)/x - & + & 0.01252108546759328*Log(rhoa) + 0.005806550506277202*t*Log(rhoa) - & + & 0.0001016735312443444*t**2*Log(rhoa) + 2.881946187214505e-7*t**3*Log(rhoa) + & + & (0.0942243379396279*Log(rhoa))/x - 0.0385459592773097*Log(rh)*Log(rhoa) - & + & 0.0006723156277391984*t*Log(rh)*Log(rhoa) + 2.602884877659698e-6*t**2*Log(rh)*Log(rhoa) + & + & 1.194163699688297e-8*t**3*Log(rh)*Log(rhoa) - (0.00851515345806281*Log(rh)*Log(rhoa))/x - & + & 0.01837488495738111*Log(rh)**2*Log(rhoa) + 0.0001720723574407498*t*Log(rh)**2*Log(rhoa) - & + & 3.717657974086814e-7*t**2*Log(rh)**2*Log(rhoa) - & + & 5.148746022615196e-10*t**3*Log(rh)**2*Log(rhoa) + & + & (0.0002686602132926594*Log(rh)**2*Log(rhoa))/x - 0.06199739728812199*Log(rhoa)**2 + & + & 0.000906958053583576*t*Log(rhoa)**2 - 9.11727926129757e-7*t**2*Log(rhoa)**2 - & + & 5.367963396508457e-9*t**3*Log(rhoa)**2 - (0.007742343393937707*Log(rhoa)**2)/x + & + & 0.0121827103101659*Log(rh)*Log(rhoa)**2 - 0.0001066499571188091*t*Log(rh)*Log(rhoa)**2 + & + & 2.534598655067518e-7*t**2*Log(rh)*Log(rhoa)**2 - & + & 3.635186504599571e-10*t**3*Log(rh)*Log(rhoa)**2 + & + & (0.0006100650851863252*Log(rh)*Log(rhoa)**2)/x + 0.0003201836700403512*Log(rhoa)**3 - & + & 0.0000174761713262546*t*Log(rhoa)**3 + 6.065037668052182e-8*t**2*Log(rhoa)**3 - & + & 1.421771723004557e-11*t**3*Log(rhoa)**3 + (0.0001357509859501723*Log(rhoa)**3)/x + ntot=exp(ntot) + + rc=exp(-1.6524245+0.42316402*x+0.33466487*log(ntot)) !nm + + IF (jnuc < 1.e-7) THEN +! print *,'Warning (ilon=',ilon,'ilev=',ilev'): nucleation rate < 1e-7/cm3s, using 0.0/cm3s,' + jnuc=0.0 + ENDIF + + IF (jnuc > 1.e10) THEN +! print *,'Warning (ilon=',ilon,'ilev=',ilev, & +! & '): nucleation rate > 1e10/cm3s, using 1e10/cm3s' + jnuc=1.e10 + ENDIF + + rhotres=exp( -279.2430007512709 + 11.73439886096903*rh + 22700.92970508331/t & + & - (1088.644983466801*rh)/t + 1.144362942094912*t & + & - 0.03023314602163684*rh*t - 0.001302541390154324*t**2 & + & - 6.386965238433532*Log(rh) + (854.980361026715*Log(rh))/t & + & + 0.00879662256826497*t*Log(rh)) !1/cm3 + + RETURN + +END SUBROUTINE binapara + +END MODULE nucleation_tstep_mod Index: LMDZ5/trunk/libf/phylmd/StratAer/ocs_to_so2.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/StratAer/ocs_to_so2.F90 (revision 2690) +++ LMDZ5/trunk/libf/phylmd/StratAer/ocs_to_so2.F90 (revision 2690) @@ -0,0 +1,44 @@ +SUBROUTINE OCS_TO_SO2(pdtphys,tr_seri,t_seri,pplay,paprs,sh,OCS_lifetime,is_strato,ocs_convert_sub) + + USE dimphy, ONLY : klon,klev + USE aerophys + USE infotrac + USE YOMCST, ONLY : RG + + IMPLICIT NONE + + !-------------------------------------------------------- + ! transfer variables when calling this routine + REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) + REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique + REAL,DIMENSION(klon,klev),INTENT(IN) :: OCS_lifetime + LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato + + ! local variables + INTEGER :: i,j,k,nb,ilon,ilev + REAL,DIMENSION(klon,klev) :: ocs_convert_sub ! OCS converted to SO2 [kg(S)/m2/layer/s] + +!--convert OCS to SO2 + DO ilon=1, klon + DO ilev=1, klev + !only in the stratosphere + IF (is_strato(ilon,ilev)) THEN + IF (OCS_lifetime(ilon,ilev).GT.0.0) THEN + ocs_convert_sub(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1-exp(-pdtphys/OCS_lifetime(ilon,ilev))) + ELSE + ocs_convert_sub(ilon,ilev)=0.0 + ENDIF + tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - ocs_convert_sub(ilon,ilev) + tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*ocs_convert_sub(ilon,ilev) + !convert ocs_convert from kg(OCS)/kgA to kg(S)/m2/layer/s for saving as diagnostic + ocs_convert_sub(ilon,ilev)=ocs_convert_sub(ilon,ilev)*mSatom/mOCSmol* & + & (paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys + ENDIF + ENDDO + ENDDO + +END SUBROUTINE OCS_TO_SO2 Index: LMDZ5/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90 (revision 2690) +++ LMDZ5/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90 (revision 2690) @@ -0,0 +1,44 @@ +SUBROUTINE SO2_TO_H2SO4(pdtphys,tr_seri,t_seri,pplay,paprs,sh,SO2_lifetime,is_strato,sulf_convert_sub) + + USE dimphy, ONLY : klon,klev + USE aerophys + USE infotrac + USE YOMCST, ONLY : RG + + IMPLICIT NONE + + !-------------------------------------------------------- + + ! transfer variables when calling this routine + REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) + REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique + REAL,DIMENSION(klon,klev),INTENT(IN) :: SO2_lifetime + LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato + + ! local variables in coagulation routine + INTEGER :: i,j,k,nb,ilon,ilev + REAL,DIMENSION(klon,klev) :: sulf_convert_sub ! SO2 converted to H2SO4 [kg(S)/m2/layer/s] + +!--convert SO2 to H2SO4 + DO ilon=1, klon + DO ilev=1, klev + !only in the stratosphere + IF (is_strato(ilon,ilev)) THEN + IF (SO2_lifetime(ilon,ilev).GT.0.0) THEN + sulf_convert_sub(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1-exp(-pdtphys/SO2_lifetime(ilon,ilev))) + ELSE + sulf_convert_sub(ilon,ilev)=0.0 + ENDIF + tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - sulf_convert_sub(ilon,ilev) + tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*sulf_convert_sub(ilon,ilev) + !convert sulf_convert from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic + sulf_convert_sub(ilon,ilev)=sulf_convert_sub(ilon,ilev)*mSatom/mSO2mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys + ENDIF + ENDDO + ENDDO + +END SUBROUTINE SO2_TO_H2SO4 Index: LMDZ5/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90 (revision 2690) +++ LMDZ5/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90 (revision 2690) @@ -0,0 +1,624 @@ +MODULE sulfate_aer_mod + +! microphysical routines based on UPMC aerosol model by Slimane Bekki +! adapted for stratospheric sulfate aerosol in LMDZ by Christoph Kleinschmitt + +CONTAINS + +!******************************************************************** + SUBROUTINE STRACOMP(sh,t_seri,pplay) + +! AEROSOL H2SO4 WEIGHT FRACTION AS A FUNCTION OF PH2O AND TEMPERATURE +! ---------------------------------------------------------------- +! INPUT: +! H2O: VMR of H2O +! t_seri: temperature (K) +! PMB: pressure (mb) +! klon: number of latitude bands in the model domain +! klev: number of altitude bands in the model domain +! for IFS: perhaps add another dimension for longitude +! +! OUTPUT: +! R2SO4: aerosol H2SO4 weight fraction (percent) + + USE dimphy, ONLY : klon,klev + USE aerophys + USE phys_local_var_mod, ONLY: R2SO4 + + IMPLICIT NONE + + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique + + REAL PMB(klon,klev), H2O(klon,klev) +! +! working variables + INTEGER I,J,K + REAL TP, PH2O, VAL, A, B +! local variables to be saved on exit + INTEGER INSTEP + INTEGER, PARAMETER :: N=16, M=28 + DATA INSTEP/0/ + REAL F(N,M) + REAL XC(N) + REAL YC(M) + REAL XC1, XC16, YC1, YC28 +! + SAVE INSTEP,F,XC,YC,XC1,XC16,YC1,YC28 + +! convert pplay (in Pa) to PMB (in mb) + PMB(:,:)=pplay(:,:)/100.0 + +! convert specific humidity sh (in kg/kg) to VMR H2O + H2O(:,:)=sh(:,:)*mAIRmol/mH2Omol + + IF(INSTEP.EQ.0) THEN + + INSTEP=1 + XC(1)=0.01 + XC(2)=0.1 + XC(3)=0.5 + XC(4)=1.0 + XC(5)=1.5 + XC(6)=2.0 + XC(7)=3.0 + XC(8)=5.0 + XC(9)=6.0 + XC(10)=8.0 + XC(11)=10.0 + XC(12)=12.0 + XC(13)=15.0 + XC(14)=20.0 + XC(15)=30.0 + XC(16)=100.0 +! + YC(1)=175.0 + DO I=2,28 + YC(I)=YC(I-1)+5.0 + ENDDO + +! CONVERSION mb IN 1.0E-4mB + DO I=1,16 + XC(I)=XC(I)*1.0E-4 + ENDDO +! + XC1=XC(1)+1.E-10 + XC16=XC(16)-1.E-8 + YC1=YC(1)+1.E-5 + YC28=YC(28)-1.E-5 + + F(6,4)=43.45 + F(6,5)=53.96 + F(6,6)=60.62 + F(6,7)=65.57 + F(6,8)=69.42 + F(6,9)=72.56 + F(6,10)=75.17 + F(6,11)=77.38 + F(6,12)=79.3 + F(6,13)=80.99 + F(6,14)=82.5 + F(6,15)=83.92 + F(6,16)=85.32 + F(6,17)=86.79 + F(6,18)=88.32 +! +! ADD FACTOR BECAUSE THE SLOP IS TOO IMPORTANT +! NOT FOR THIS ONE BUT THE REST +! LOG DOESN'T WORK + A=(F(6,5)-F(6,4))/( (YC(5)-YC(4))*2.0) + B=-A*YC(4) + F(6,4) + F(6,1)=A*YC(1) + B + F(6,2)=A*YC(2) + B + F(6,3)=A*YC(3) + B +! + F(7,4)=37.02 + F(7,5)=49.46 + F(7,6)=57.51 + F(7,7)=63.12 + F(7,8)=67.42 + F(7,9)=70.85 + F(7,10)=73.70 + F(7,11)=76.09 + F(7,12)=78.15 + F(7,13)=79.96 + F(7,14)=81.56 + F(7,15)=83.02 + F(7,16)=84.43 + F(7,17)=85.85 + F(7,18)=87.33 +! + A=(F(7,5)-F(7,4))/( (YC(5)-YC(4))*2.0) + B=-A*YC(4) + F(7,4) + F(7,1)=A*YC(1) + B + F(7,2)=A*YC(2) + B + F(7,3)=A*YC(3) + B +! + F(8,4)=25.85 + F(8,5)=42.26 + F(8,6)=52.78 + F(8,7)=59.55 + F(8,8)=64.55 + F(8,9)=68.45 + F(8,10)=71.63 + F(8,11)=74.29 + F(8,12)=76.56 + F(8,13)=78.53 + F(8,14)=80.27 + F(8,15)=81.83 + F(8,16)=83.27 + F(8,17)=84.67 + F(8,18)=86.10 +! + A=(F(8,5)-F(8,4))/( (YC(5)-YC(4))*2.5 ) + B=-A*YC(4) + F(8,4) + F(8,1)=A*YC(1) + B + F(8,2)=A*YC(2) + B + F(8,3)=A*YC(3) + B +! + F(9,4)=15.38 + F(9,5)=39.35 + F(9,6)=50.73 + F(9,7)=58.11 + F(9,8)=63.41 + F(9,9)=67.52 + F(9,10)=70.83 + F(9,11)=73.6 + F(9,12)=75.95 + F(9,13)=77.98 + F(9,14)=79.77 + F(9,15)=81.38 + F(9,16)=82.84 + F(9,17)=84.25 + F(9,18)=85.66 +! + A=(F(9,5)-F(9,4))/( (YC(5)-YC(4))*7.0) + B=-A*YC(4) + F(9,4) + F(9,1)=A*YC(1) + B + F(9,2)=A*YC(2) + B + F(9,3)=A*YC(3) + B +! + F(10,4)=0.0 + F(10,5)=34.02 + F(10,6)=46.93 + F(10,7)=55.61 + F(10,8)=61.47 + F(10,9)=65.94 + F(10,10)=69.49 + F(10,11)=72.44 + F(10,12)=74.93 + F(10,13)=77.08 + F(10,14)=78.96 + F(10,15)=80.63 + F(10,16)=82.15 + F(10,17)=83.57 + F(10,18)=84.97 +! + A=(F(10,6)-F(10,5))/( (YC(6)-YC(5))*1.5) + B=-A*YC(5) + F(10,5) + F(10,1)=A*YC(1) + B + F(10,2)=A*YC(2) + B + F(10,3)=A*YC(3) + B + F(10,4)=A*YC(4) + B +! + F(11,4)=0.0 + F(11,5)=29.02 + F(11,6)=43.69 + F(11,7)=53.44 + F(11,8)=59.83 + F(11,9)=64.62 + F(11,10)=68.39 + F(11,11)=71.48 + F(11,12)=74.10 + F(11,13)=76.33 + F(11,14)=78.29 + F(11,15)=80.02 + F(11,16)=81.58 + F(11,17)=83.03 + F(11,18)=84.44 +! + A=(F(11,6)-F(11,5))/( (YC(6)-YC(5))*2.5 ) + B=-A*YC(5) + F(11,5) + F(11,1)=A*YC(1) + B + F(11,2)=A*YC(2) + B + F(11,3)=A*YC(3) + B + F(11,4)=A*YC(4) + B +! + F(12,4)=0.0 + F(12,5)=23.13 + F(12,6)=40.86 + F(12,7)=51.44 + F(12,8)=58.38 + F(12,9)=63.47 + F(12,10)=67.43 + F(12,11)=70.66 + F(12,12)=73.38 + F(12,13)=75.70 + F(12,14)=77.72 + F(12,15)=79.51 + F(12,16)=81.11 + F(12,17)=82.58 + F(12,18)=83.99 +! + A=(F(12,6)-F(12,5))/( (YC(6)-YC(5))*3.5 ) + B=-A*YC(5) + F(12,5) + F(12,1)=A*YC(1) + B + F(12,2)=A*YC(2) + B + F(12,3)=A*YC(3) + B + F(12,4)=A*YC(4) + B +! + F(13,4)=0.0 + F(13,5)=0.0 + F(13,6)=36.89 + F(13,7)=48.63 + F(13,8)=56.46 + F(13,9)=61.96 + F(13,10)=66.19 + F(13,11)=69.6 + F(13,12)=72.45 + F(13,13)=74.89 + F(13,14)=76.99 + F(13,15)=78.85 + F(13,16)=80.50 + F(13,17)=82.02 + F(13,18)=83.44 +! + A=(F(13,7)-F(13,6))/( (YC(7)-YC(6))*2.0) + B=-A*YC(6) + F(13,6) + F(13,1)=A*YC(1) + B + F(13,2)=A*YC(2) + B + F(13,3)=A*YC(3) + B + F(13,4)=A*YC(4) + B + F(13,5)=A*YC(5) + B +! + F(14,4)=0.0 + F(14,5)=0.0 + F(14,6)=30.82 + F(14,7)=44.49 + F(14,8)=53.69 + F(14,9)=59.83 + F(14,10)=64.47 + F(14,11)=68.15 + F(14,12)=71.19 + F(14,13)=73.77 + F(14,14)=76.0 + F(14,15)=77.95 + F(14,16)=79.69 + F(14,17)=81.26 + F(14,18)=82.72 +! + A=(F(14,7)-F(14,6))/( (YC(7)-YC(6))*2.5 ) + B=-A*YC(6) + F(14,6) + F(14,1)=A*YC(1) + B + F(14,2)=A*YC(2) + B + F(14,3)=A*YC(3) + B + F(14,4)=A*YC(4) + B + F(14,5)=A*YC(5) + B +! + F(15,4)=0.0 + F(15,5)=0.0 + F(15,6)=0.0 + F(15,7)=37.71 + F(15,8)=48.49 + F(15,9)=56.40 + F(15,10)=61.75 + F(15,11)=65.89 + F(15,12)=69.25 + F(15,13)=72.07 + F(15,14)=74.49 + F(15,15)=76.59 + F(15,16)=78.45 + F(15,17)=80.12 + F(15,18)=81.64 +! + A=(F(15,8)-F(15,7))/( (YC(8)-YC(7))*1.5) + B=-A*YC(7) + F(15,7) + F(15,1)=A*YC(1) + B + F(15,2)=A*YC(2) + B + F(15,3)=A*YC(3) + B + F(15,4)=A*YC(4) + B + F(15,5)=A*YC(5) + B + F(15,6)=A*YC(6) + B + +! SUPPOSE THAT AT GIVEN AND PH2O<2mB, +! %H2SO4 = A *LOG(PH2O) +B +! XC(1-5) :EXTENSION LEFT (LOW H2O) + DO J=1,18 + A=(F(6,J)-F(7,J))/(LOG(XC(6))-LOG(XC(7))) + B=-A*LOG(XC(6)) + F(6,J) + DO K=1,5 + F(K,J)=A*LOG(XC(K)) + B + ENDDO + ENDDO + +! XC(16) :EXTENSION RIGHT (HIGH H2O) + DO J=1,18 + A=(F(15,J)-F(14,J))/(XC(15)-XC(14)) + B=-A*XC(15) + F(15,J) + F(16,J)=A*XC(16) + B +! F(16,2)=1.0 + ENDDO + +! YC(16-25) :EXTENSION DOWN (HIGH T) + DO I=1,16 + A=(F(I,18)-F(I,17))/(YC(18)-YC(17)) + B=-A*YC(18) + F(I,18) + DO K=19,28 + F(I,K)=A*YC(K) + B + ENDDO + ENDDO + +! MANUAL CORRECTIONS + DO J=1,10 + F(1,J)=94.0 + ENDDO + + DO J=1,6 + F(2,J)=77.0 +REAL(J) + ENDDO + + DO J=1,7 + F(16,J)=9.0 + ENDDO + + DO I=1,16 + DO J=1,28 + IF (F(I,J).LT.9.0) F(I,J)=30.0 + IF (F(I,J).GT.99.99) F(I,J)=99.99 + ENDDO + ENDDO + + ENDIF + + DO I=1,klon + DO J=1,klev + TP=t_seri(I,J) + IF (TP.LT.175.1) TP=175.1 +! Partial pressure of H2O (mb) + PH2O =PMB(I,J)*H2O(I,J) + IF (PH2O.LT.XC1) THEN + R2SO4(I,J)=99.99 +! PH2O=XC(1)+1.0E-10 + ELSE + IF (PH2O.GT.XC16) PH2O=XC16 +! SIMPLE LINEAR INTERPOLATIONS + CALL FIND(PH2O,TP,XC,YC,F,VAL,N,M) + IF (PMB(I,J).GE.10.0.AND.VAL.LT.60.0) VAL=60.0 + R2SO4(I,J)=VAL + ENDIF + ENDDO + ENDDO + + END SUBROUTINE + +!**************************************************************** + SUBROUTINE STRAACT(ACTSO4) + +! H2SO4 ACTIVITY (GIAUQUE) AS A FUNCTION OF H2SO4 WP +! ---------------------------------------- +! INPUT: +! H2SO4: VMR of H2SO4 +! klon: number of latitude bands in the model domain +! klev: number of altitude bands in the model domain +! for IFS: perhaps add another dimension for longitude +! +! OUTPUT: +! ACTSO4: H2SO4 activity (percent) + + USE dimphy, ONLY : klon,klev + USE phys_local_var_mod, ONLY: R2SO4 + + IMPLICIT NONE + + REAL ACTSO4(klon,klev) + +! Working variables + INTEGER NN,I,J,JX,JX1 + REAL TC,TB,TA,XT + PARAMETER (NN=109) + REAL XC(NN), X(NN) + +! H2SO4 activity + DATA X/ & + & 0.0,0.25,0.78,1.437,2.19,3.07,4.03,5.04,6.08 & + & ,7.13,8.18,14.33,18.59,28.59,39.17,49.49 & + & ,102.4,157.8,215.7,276.9,341.6,409.8,481.5,556.6 & + & ,635.5,719.,808.,902.,1000.,1103.,1211.,1322.,1437.,1555. & + & ,1677.,1800.,1926.,2054.,2183.,2312.,2442.,2572.,2701.,2829. & + & ,2955.,3080.,3203.,3325.,3446.,3564.,3681.,3796.,3910.,4022. & + & ,4134.,4351.,4564.,4771.,4974.,5171.,5364.,5551.,5732.,5908. & + & ,6079.,6244.,6404.,6559.,6709.,6854.,6994.,7131.,7264.,7393. & + & ,7520.,7821.,8105.,8373.,8627.,8867.,9093.,9308.,9511.,9703. & + & ,9885.,10060.,10225.,10535.,10819.,11079.,11318.,11537. & + & ,11740.,12097.,12407.,12676.,12915.,13126.,13564.,13910. & + & ,14191.,14423.,14617.,14786.,10568.,15299.,15491.,15654. & + & ,15811./ +! H2SO4 weight fraction (percent) + DATA XC/ & + & 100.0,99.982,99.963,99.945,99.927,99.908,99.890,99.872 & + & ,99.853,99.835,99.817,99.725,99.634,99.452,99.270 & + & ,99.090,98.196,97.319,96.457,95.610,94.777,93.959,93.156 & + & ,92.365,91.588,90.824,90.073,89.334,88.607,87.892,87.188 & + & ,86.495,85.814,85.143,84.482,83.832,83.191,82.560,81.939 & + & ,81.327,80.724,80.130,79.545,78.968,78.399,77.839,77.286 & + & ,76.741,76.204,75.675,75.152,74.637,74.129,73.628,73.133 & + & ,72.164,71.220,70.300,69.404,68.530,67.678,66.847,66.037 & + & ,65.245,64.472,63.718,62.981,62.261,61.557,60.868,60.195 & + & ,59.537,58.893,58.263,57.646,56.159,54.747,53.405,52.126 & + & ,50.908,49.745,48.634,47.572,46.555,45.580,44.646,43.749 & + & ,42.059,40.495,39.043,37.691,36.430,35.251,33.107,31.209 & + & ,29.517,27.999,26.629,23.728,21.397,19.482,17.882,16.525 & + & ,15.360,13.461,11.980,10.792,9.819,8.932/ + + DO I=1,klon + DO J=1,klev +! HERE LINEAR INTERPOLATIONS + XT=R2SO4(I,J) + CALL POSACT(XT,XC,NN,JX) + JX1=JX+1 + IF(JX.EQ.0) THEN + ACTSO4(I,J)=0.0 + ELSE IF(JX.GE.NN) THEN + ACTSO4(I,J)=15811.0 + ELSE + TC=XT -XC(JX) + TB=X(JX1) -X(JX) + TA=XC(JX1) -XC(JX) + TA=TB/TA + ACTSO4(I,J)=X(JX) + TA*TC + ENDIF + ENDDO + ENDDO + + END SUBROUTINE + +!**************************************************************** + SUBROUTINE DENH2SA(t_seri) + +! AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T +! --------------------------------------------- +! VERY ROUGH APPROXIMATION (SEE FOR WATER IN HANDBOOK +! LINEAR 2% FOR 30 DEGREES with RESPECT TO WATER) +! +! INPUT: +! R2SO4: aerosol H2SO4 weight fraction (percent) +! t_seri: temperature (K) +! klon: number of latitude bands in the model domain +! klev: number of altitude bands in the model domain +! for IFS: perhaps add another dimension for longitude +! +! OUTPUT: +! DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume) +! + USE dimphy, ONLY : klon,klev + USE phys_local_var_mod, ONLY: R2SO4, DENSO4 + + IMPLICIT NONE + + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + + INTEGER I,J + +! Loop on model domain (2 dimension for UPMC model; 3 for IFS) + DO I=1,klon + DO J=1,klev +! RO AT 20C + DENSO4(I,J)=0.78681252E-5*R2SO4(I,J)*R2SO4(I,J)+ 0.82185978E-2*R2SO4(I,J)+0.97968381 + DENSO4(I,J)=DENSO4(I,J)* ( 1.0 - (t_seri(I,J)-293.0)*0.02/30.0 ) + ENDDO + ENDDO + + END SUBROUTINE + +!*********************************************************** + SUBROUTINE FIND(X,Y,XC,YC,F,VAL,N,M) +! +! BI-LINEAR INTERPOLATION + +! INPUT: +! X: Partial pressure of H2O (mb) +! Y: temperature (K) +! XC: Table partial pressure of H2O (mb) +! YC: Table temperature (K) +! F: Table aerosol H2SO4 weight fraction=f(XC,YC) (percent) +! +! OUTPUT: +! VAL: aerosol H2SO4 weight fraction (percent) + + IMPLICIT NONE + + INTEGER N,M + REAL X,Y,XC(N),YC(M),F(N,M),VAL +! +! working variables + INTEGER IERX,IERY,JX,JY,JXP1,JYP1 + REAL SXY,SX1Y,SX1Y1,SXY1,TA,TB,T,UA,UB,U + + IERX=0 + IERY=0 + CALL POSITION(XC,X,N,JX,IERX) + CALL POSITION(YC,Y,M,JY,IERY) + + IF(JX.EQ.0.OR.IERY.EQ.1) THEN + VAL=99.99 + RETURN + ENDIF + + IF(JY.EQ.0.OR.IERX.EQ.1) THEN + VAL=9.0 + RETURN + ENDIF + + JXP1=JX+1 + JYP1=JY+1 + SXY=F(JX, JY ) + SX1Y=F(JXP1,JY ) + SX1Y1=F(JXP1,JYP1) + SXY1=F(JX, JYP1) + +! x-slope. + TA=X -XC(JX) + TB=XC(JXP1)-XC(JX) + T=TA/TB + +! y-slope. + UA=Y -YC(JY) + UB=YC(JYP1)-YC(JY) + U=UA/UB + +! Use bilinear interpolation to determine function at point X,Y. + VAL=(1.-T)*(1.-U)*SXY + T*(1.0-U)*SX1Y + T*U*SX1Y1 + (1.0-T)*U*SXY1 + + IF(VAL.LT.9.0) VAL=9.0 + IF(VAL.GT.99.99) VAL=99.99 + + RETURN + END SUBROUTINE +!**************************************************************** + SUBROUTINE POSITION(XC,X,N,JX,IER) + + IMPLICIT NONE + + INTEGER N,JX,IER,I + REAL X,XC(N) + + IER=0 + IF(X.LT.XC(1)) THEN + JX=0 + ELSE + DO 10 I=1,N + IF (X.LT.XC(I)) GO TO 20 + 10 CONTINUE + IER=1 + 20 JX=I-1 + ENDIF + + RETURN + END SUBROUTINE +!******************************************************************** + SUBROUTINE POSACT(XT,X,N,JX) + +! POSITION OF XT IN THE ARRAY X +! ----------------------------------------------- + + IMPLICIT NONE + + INTEGER N + REAL XT,X(N) +! Working variables + INTEGER JX,I + + IF(XT.GT.X(1)) THEN + JX=0 + ELSE + DO 10 I=1,N + IF (XT.GT.X(I)) GO TO 20 + 10 CONTINUE + 20 JX=I + ENDIF + + RETURN + END SUBROUTINE + +END MODULE sulfate_aer_mod Index: LMDZ5/trunk/libf/phylmd/StratAer/traccoag_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/StratAer/traccoag_mod.F90 (revision 2690) +++ LMDZ5/trunk/libf/phylmd/StratAer/traccoag_mod.F90 (revision 2690) @@ -0,0 +1,321 @@ +MODULE traccoag_mod +! +! This module calculates the concentration of aerosol particles in certain size bins +! considering coagulation and sedimentation. +! +CONTAINS + + SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, & + presnivs, xlat, xlon, pphis, pphi, & + t_seri, pplay, paprs, sh, rh , & + tr_seri) + + USE phys_local_var_mod, ONLY: mdw, sulf_convert, sulf_nucl, sulf_cond_evap, & + & sfluxaer, ocs_convert, R2SO4, DENSO4, f_r_wet, SO2_backgr_tend, OCS_backgr_tend, & + & OCS_lifetime, SO2_lifetime, surf_PM25_sulf + + USE dimphy + USE infotrac + USE aerophys + USE geometry_mod, ONLY : cell_area + USE mod_grid_phy_lmdz + USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root + USE mod_phys_lmdz_para, only: gather, scatter + USE phys_cal_mod + USE sulfate_aer_mod + USE nucleate_mod + USE phys_local_var_mod, ONLY: stratomask + USE sulfate_aer_mod, ONLY : stracomp, denh2sa + USE YOMCST + + IMPLICIT NONE + +! Input argument +!--------------- + REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) + REAL,INTENT(IN) :: gmtime ! Heure courante + LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique + INTEGER,INTENT(IN) :: julien ! Jour julien + + REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA) + REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point + REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point + REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol + REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche + + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique + REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative + +! Output argument +!---------------- + REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] + +! Local variables +!---------------- +! flag for sulfur emission scenario: (0) background aerosol ; (1) volcanic eruption ; (2) stratospheric aerosol injections (SAI) + INTEGER,PARAMETER :: flag_sulf_emit=2 +! +!--flag_sulf_emit=1 --example Pinatubo + INTEGER,PARAMETER :: year_emit_vol=1991 ! year of emission date + INTEGER,PARAMETER :: mth_emit_vol=6 ! month of emission date + INTEGER,PARAMETER :: day_emit_vol=15 ! day of emission date + REAL,PARAMETER :: m_aer_emiss_vol=7.e9 ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2) + REAL,PARAMETER :: altemiss_vol=17.e3 ! emission altitude in m + REAL,PARAMETER :: sigma_alt_vol=1.e3 ! standard deviation of emission altitude in m + REAL,PARAMETER :: xlat_vol=15.14 ! latitude of volcano in degree + REAL,PARAMETER :: xlon_vol=120.35 ! longitude of volcano in degree + +!--flag_sulf_emit=2 --SAI + REAL,PARAMETER :: m_aer_emiss_sai=1.e10 ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS + REAL,PARAMETER :: altemiss_sai=17.e3 ! emission altitude in m + REAL,PARAMETER :: sigma_alt_sai=1.e3 ! standard deviation of emission altitude in m + REAL,PARAMETER :: xlat_sai=0.0 ! latitude of SAI in degree + REAL,PARAMETER :: xlon_sai=120.35 ! longitude of SAI in degree + +!--be careful - this needs to be changed with resolution - here for 96x95 + REAL,PARAMETER :: dlat=0.9474 ! d latitude in degree + REAL,PARAMETER :: dlon=1.875 ! d longitude in degree + +!--other local variables + INTEGER :: it, k, i, ilon, ilev, itime, i_int + LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below + REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg] + REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA] + REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m + REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer + REAL :: f_lay_sum ! sum of layer emission fractions + REAL :: alt_lower ! integral boundary for Gaussian emission profile in altitude + REAL :: alt_upper ! integral boundary for Gaussian emission profile in altitude + INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile + REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m] + REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m] + REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m] + REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg] + REAL :: zrho ! Density of air [kg/m3] + REAL :: zdz ! thickness of atm. model layer in m + REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction + REAL,DIMENSION(klon,klev) :: sulf_convert_sub ! SO2 converted to H2SO4 [kg(S)/m2/layer/s] + REAL,DIMENSION(klon,klev) :: sulf_nucl_sub ! H2SO4 converted from gas to particles [kg(S)/m2/layer/s] + REAL,DIMENSION(klon,klev) :: ocs_convert_sub ! OCS converted to SO2 [kg(S)/m2/layer/s] + REAL,DIMENSION(klon,klev) :: SO2_backgr_tend_sub ! SO2 from background [kg(S)/m2/layer/s] + REAL,DIMENSION(klon,klev) :: OCS_backgr_tend_sub ! OCS from background [kg(S)/m2/layer/s] + + CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'deb SO2') + DO it=1, nbtr_bin + CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'deb bin ') + ENDDO + + IF (is_mpi_root) THEN + PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour + ENDIF + + DO it=1, nbtr_bin + r_bin(it)=mdw(it)/2. + ENDDO + +!--set boundaries of size bins + DO it=1, nbtr_bin + IF (it.EQ.1) THEN + r_upper(it)=sqrt(r_bin(it+1)*r_bin(it)) + r_lower(it)=r_bin(it)**2./r_upper(it) + ELSEIF (it.EQ.nbtr_bin) THEN + r_lower(it)=sqrt(r_bin(it)*r_bin(it-1)) + r_upper(it)=r_bin(it)**2./r_lower(it) + ELSE + r_lower(it)=sqrt(r_bin(it)*r_bin(it-1)) + r_upper(it)=sqrt(r_bin(it+1)*r_bin(it)) + ENDIF + ENDDO + + IF (debutphy .and. is_mpi_root) THEN + DO it=1, nbtr_bin + PRINT *,'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' + ENDDO + ENDIF + +!--initialising logical is_strato from stratomask + is_strato(:,:)=.FALSE. + WHERE (stratomask(:,:).GT.0.5) is_strato(:,:)=.TRUE. + +! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) +! H2SO4 mass fraction in aerosol (%) + CALL stracomp(sh,t_seri,pplay) + +! aerosol density (gr/cm3) + CALL denh2sa(t_seri) + +! compute factor for converting dry to wet radius (for every grid box) + f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.) + +!--calculate mass of air in every grid box + DO ilon=1, klon + DO ilev=1, klev + m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1)) / RG * cell_area(ilon) + ENDDO + ENDDO + + IF (debutphy) THEN + CALL gather(tr_seri, tr_seri_glo) + IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN +!--initialising tracer concentrations to zero + DO it=1, nbtr + tr_seri(:,:,it)=0.0 + ENDDO + ENDIF + ENDIF + +!--sulfur emission, depending on chosen scenario (flag_sulf_emit) + SELECT CASE(flag_sulf_emit) + + CASE(0) ! background aerosol + ! do nothing (no emission) + + CASE(1) ! volcanic eruption + !--only emit on day of eruption + ! stretch emission over one day of Pinatubo eruption + IF (year_cur==year_emit_vol.AND.mth_cur==mth_emit_vol.AND.day_cur==day_emit_vol) THEN +! + DO i=1,klon + !Pinatubo eruption at 15.14N, 120.35E + IF ( xlat(i).GT.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. & + xlon(i).GT.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN +! compute altLMDz + altLMDz(:)=0.0 + DO k=1, klev + zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 + zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG !thickness of layer in m + altLMDz(k+1)=altLMDz(k)+zdz + ENDDO + !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) + f_lay_sum=0.0 + DO k=1, klev + f_lay_emiss(k)=0.0 + alt_lower=altLMDz(k) + DO i_int=1, n_int_alt + alt_upper=alt_lower+(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) + f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol)* & + & (exp(-0.5*((altemiss_vol-alt_lower)/sigma_alt_vol)**2.)+ & + & exp(-0.5*((alt_upper-altemiss_vol)/sigma_alt_vol)**2.))* & + & (alt_upper-alt_lower)/2. + ENDDO + f_lay_sum=f_lay_sum+f_lay_emiss(k) + ENDDO + !correct for step integration error + f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum + !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) + !vertically distributed emission + DO k=1, klev + tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+ & + & m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k) & + & /(1.*86400./pdtphys) ! stretch emission over one day of Pinatubo eruption + ENDDO + ENDIF ! emission grid cell + ENDDO ! klon loop + ENDIF ! emission period + + CASE(2) ! stratospheric aerosol injections (SAI) +! + DO i=1,klon +! SAI standard scenario with continuous emission from 1 grid point at the equator +! SAI emission on single month +! IF ((mth_cur==4 .AND. & +! SAI continuous emission o + IF ( xlat(i).GT.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. & + & xlon(i).GT.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN +! compute altLMDz + altLMDz(:)=0.0 + DO k=1, klev + zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 + zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG !thickness of layer in m + altLMDz(k+1)=altLMDz(k)+zdz + ENDDO + !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) + f_lay_sum=0.0 + DO k=1, klev + f_lay_emiss(k)=0.0 + alt_lower=altLMDz(k) + DO i_int=1, n_int_alt + alt_upper=alt_lower+(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) + f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & + & (exp(-0.5*((altemiss_sai-alt_lower)/sigma_alt_sai)**2.)+ & + & exp(-0.5*((alt_upper-altemiss_sai)/sigma_alt_sai)**2.))* & + & (alt_upper-alt_lower)/2. + ENDDO + f_lay_sum=f_lay_sum+f_lay_emiss(k) + ENDDO + !correct for step integration error + f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum + !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) + !vertically distributed emission + DO k=1, klev + tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+ & + & m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k) & + & /(360.*86400./pdtphys) ! stretch emission over whole year (360d) +! & /(60.*86400./pdtphys) ! stretch emission over 2 months (seasonal emission) +! & /7. ! distribute equally over 7 emission grid points + ENDDO +! !emission as monodisperse particles with 0.1um dry radius (BIN21) +! !vertically distributed emission +! DO k=1, klev +! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+ & +! & m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k) & +! & /(360.*86400./pdtphys) & ! stretch emission over whole year (360d) +! & /7. ! distribute equally over 7 emission grid points +! ENDDO + ENDIF ! emission grid cell + ENDDO ! klon loop + + END SELECT ! emission scenario (flag_sulf_emit) + +!--read background concentrations of OCS and SO2 and lifetimes from input file + SO2_backgr_tend_sub(:,:)=0.0 + OCS_backgr_tend_sub(:,:)=0.0 + CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri,SO2_backgr_tend_sub, & + & OCS_backgr_tend_sub,SO2_lifetime,OCS_lifetime) + SO2_backgr_tend(:,:)=SO2_backgr_tend_sub(:,:) + OCS_backgr_tend(:,:)=OCS_backgr_tend_sub(:,:) + +!--convert OCS to SO2 in the stratosphere + ocs_convert_sub(:,:)=0.0 + CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,sh,OCS_lifetime,is_strato,ocs_convert_sub) + ocs_convert(:,:)=ocs_convert_sub(:,:) + +!--convert SO2 to H2SO4 + sulf_convert_sub(:,:)=0.0 + CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,sh,SO2_lifetime,is_strato,sulf_convert_sub) + sulf_convert(:,:)=sulf_convert_sub(:,:) + +!--common routine for nucleation and condensation/evaporation with adaptive timestep + CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato) + +!--call coagulation routine + CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) + +!--call sedimentation routine + CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) + +!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies + surf_PM25_sulf(:)=0.0 + DO i=1,klon + DO it=1, nbtr_bin + IF (mdw(it) .LT. 2.5e-6) THEN + !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & + !assume that particles consist of ammonium sulfate at the surface (132g/mol) and are dry at T = 20 deg. C and 50 perc. humidity + surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) & + & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 & + & *pplay(i,1)/t_seri(i,1)/RD*1e9 + ENDIF + ENDDO + ENDDO + + CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2') + DO it=1, nbtr_bin + CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ') + ENDDO + + END SUBROUTINE traccoag + +END MODULE traccoag_mod Index: LMDZ5/trunk/libf/phylmd/infotrac_phy.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/infotrac_phy.F90 (revision 2689) +++ LMDZ5/trunk/libf/phylmd/infotrac_phy.F90 (revision 2690) @@ -95,7 +95,9 @@ indnum_fn_num_,index_trac_,& niso_,ntraceurs_zone_,ntraciso_) - ! transfer information on tracers from dynamics to physics - USE print_control_mod, ONLY: prt_level, lunout - IMPLICIT NONE + + ! transfer information on tracers from dynamics to physics + USE print_control_mod, ONLY: prt_level, lunout + IMPLICIT NONE + INTEGER,INTENT(IN) :: nqtot_ INTEGER,INTENT(IN) :: nqo_ Index: LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 2689) +++ LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 2690) @@ -29,8 +29,6 @@ REAL, SAVE, ALLOCATABLE :: d_u_dyn(:,:), d_v_dyn(:,:) !$OMP THREADPRIVATE(d_u_dyn, d_v_dyn) -!!!! REAL, SAVE, ALLOCATABLE :: d_tr_dyn(:,:,:) !$OMP THREADPRIVATE(d_tr_dyn) -!!!! REAL, SAVE, ALLOCATABLE :: d_t_con(:,:),d_q_con(:,:) !$OMP THREADPRIVATE(d_t_con,d_q_con) @@ -420,4 +418,60 @@ !$OMP THREADPRIVATE(sissnow,runoff,albsol3_lic) +#ifdef CPP_StratAer +! variables for strat. aerosol CK + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: R2SO4 +!$OMP THREADPRIVATE(R2SO4) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: DENSO4 +!$OMP THREADPRIVATE(DENSO4) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: f_r_wet +!$OMP THREADPRIVATE(f_r_wet) + REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: sfluxaer +!$OMP THREADPRIVATE(sfluxaer) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: decfluxaer +!$OMP THREADPRIVATE(decfluxaer) + REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: mdw +!$OMP THREADPRIVATE(mdw) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sulf_convert +!$OMP THREADPRIVATE(sulf_convert) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sulf_nucl +!$OMP THREADPRIVATE(sulf_nucl) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sulf_cond_evap +!$OMP THREADPRIVATE(sulf_cond_evap) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ocs_convert +!$OMP THREADPRIVATE(ocs_convert) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: SO2_backgr_tend +!$OMP THREADPRIVATE(SO2_backgr_tend) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: OCS_backgr_tend +!$OMP THREADPRIVATE(OCS_backgr_tend) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: OCS_lifetime +!$OMP THREADPRIVATE(OCS_lifetime) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: SO2_lifetime +!$OMP THREADPRIVATE(SO2_lifetime) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: alpha_bin +!$OMP THREADPRIVATE(alpha_bin) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: piz_bin +!$OMP THREADPRIVATE(piz_bin) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: cg_bin +!$OMP THREADPRIVATE(cg_bin) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tau_strat_550 +!$OMP THREADPRIVATE(tau_strat_550) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tau_strat_550_lay +!$OMP THREADPRIVATE(tau_strat_550_lay) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tau_strat_1020 +!$OMP THREADPRIVATE(tau_strat_1020) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tausum_strat +!$OMP THREADPRIVATE(tausum_strat) + REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: sulf_dep_dry +!$OMP THREADPRIVATE(sulf_dep_dry) + REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: sulf_dep_wet +!$OMP THREADPRIVATE(sulf_dep_wet) + REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: surf_PM25_sulf +!$OMP THREADPRIVATE(surf_PM25_sulf) + REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: p_tropopause +!$OMP THREADPRIVATE(p_tropopause) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vsed_aer +!$OMP THREADPRIVATE(vsed_aer) +#endif + CONTAINS @@ -659,5 +713,32 @@ ALLOCATE (sissnow(klon),runoff(klon),albsol3_lic(klon)) - +#ifdef CPP_StratAer + ALLOCATE (R2SO4(klon,klev)) + ALLOCATE (DENSO4(klon,klev)) + ALLOCATE (f_r_wet(klon,klev)) + ALLOCATE (sfluxaer(klon)) + ALLOCATE (decfluxaer(klon,nbtr)) + ALLOCATE (mdw(nbtr)) + ALLOCATE (sulf_convert(klon,klev)) + ALLOCATE (sulf_nucl(klon,klev)) + ALLOCATE (sulf_cond_evap(klon,klev)) + ALLOCATE (ocs_convert(klon,klev)) + ALLOCATE (SO2_backgr_tend(klon,klev)) + ALLOCATE (OCS_backgr_tend(klon,klev)) + ALLOCATE (OCS_lifetime(klon,klev)) + ALLOCATE (SO2_lifetime(klon,klev)) + ALLOCATE (alpha_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave+nwave_lw,nbtr)) + ALLOCATE (piz_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave+nwave_lw,nbtr)) + ALLOCATE (cg_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave+nwave_lw,nbtr)) + ALLOCATE (tau_strat_550(klon,klev)) + ALLOCATE (tau_strat_550_lay(klon,klev)) + ALLOCATE (tau_strat_1020(klon,klev)) + ALLOCATE (tausum_strat(klon,3)) + ALLOCATE (sulf_dep_dry(klon)) + ALLOCATE (sulf_dep_wet(klon)) + ALLOCATE (surf_PM25_sulf(klon)) + ALLOCATE (p_tropopause(klon)) + ALLOCATE (vsed_aer(klon,klev)) +#endif END SUBROUTINE phys_local_var_init @@ -876,4 +957,34 @@ DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic) +#ifdef CPP_StratAer +! variables for strat. aerosol CK + DEALLOCATE (R2SO4) + DEALLOCATE (DENSO4) + DEALLOCATE (f_r_wet) + DEALLOCATE (sfluxaer) + DEALLOCATE (decfluxaer) + DEALLOCATE (mdw) + DEALLOCATE (sulf_convert) + DEALLOCATE (sulf_nucl) + DEALLOCATE (sulf_cond_evap) + DEALLOCATE (ocs_convert) + DEALLOCATE (SO2_backgr_tend) + DEALLOCATE (OCS_backgr_tend) + DEALLOCATE (SO2_lifetime) + DEALLOCATE (OCS_lifetime) + DEALLOCATE (alpha_bin) + DEALLOCATE (piz_bin) + DEALLOCATE (cg_bin) + DEALLOCATE (tau_strat_550) + DEALLOCATE (tau_strat_550_lay) + DEALLOCATE (tau_strat_1020) + DEALLOCATE (tausum_strat) + DEALLOCATE (sulf_dep_dry) + DEALLOCATE (sulf_dep_wet) + DEALLOCATE (surf_PM25_sulf) + DEALLOCATE (p_tropopause) + DEALLOCATE (vsed_aer) +#endif + END SUBROUTINE phys_local_var_end Index: LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 2689) +++ LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 2690) @@ -1174,4 +1174,46 @@ 'lcc', 'Cloud liquid fraction at top of cloud', '1', (/ ('', i=1, 9) /)) +#ifdef CPP_StratAer + TYPE(ctrl_out), SAVE :: o_ext_strat_550 = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'ext_strat_550', 'Strat. aerosol extinction coefficient at 550 nm', '1/m', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_ext_strat_1020 = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'ext_strat_1020', 'Strat. aerosol extinction coefficient at 1020 nm', '1/m', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_tau_strat_550 = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'OD550_strat_only', 'Stratospheric Aerosol Optical depth at 550 nm ', '1', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_tau_strat_1020 = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'OD1020_strat_only', 'Stratospheric Aerosol Optical depth at 1020 nm ', '1', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_sulf_convert = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'sulf_convert', 'SO2 mass flux converted to H2SO4', 'kg(S)/m2/layer/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_sulf_nucl = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'sulf_nucl', 'H2SO4 nucleation mass flux', 'kg(S)/m2/layer/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_sulf_cond_evap = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'sulf_cond_evap', 'H2SO4 condensation/evaporation mass flux', 'kg(S)/m2/layer/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_ocs_convert = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'ocs_convert', 'OCS mass flux converted to SO2', 'kg(S)/m2/layer/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_R2SO4 = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'R2SO4', 'H2SO4 mass fraction in aerosol', '%', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_OCS_lifetime = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'OCS_lifetime', 'OCS lifetime', 's', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_SO2_lifetime = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'SO2_lifetime', 'SO2 lifetime', 's', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_SO2_backgr_tend = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'SO2_backgr_tend', 'SO2 background tendency', 'kg(S)/m2/layer/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_OCS_backgr_tend = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'OCS_backgr_tend', 'OCS background tendency', 'kg(S)/m2/layer/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_vsed_aer = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'vsed_aer', 'Strat. aerosol sedimentation velocity (mass-weighted)', 'm/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_f_r_wet = ctrl_out((/ 1, 6, 7, 10, 10, 10, 11, 11, 11 /), & + 'f_r_wet', 'Conversion factor dry to wet aerosol radius', '-', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_sulf_dep_dry = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11 /), & + 'sulf_dep_dry', 'Sulfur dry deposition flux', 'kg(S)/m2/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_sulf_dep_wet = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11 /), & + 'sulf_dep_wet', 'Sulfur wet deposition flux', 'kg(S)/m2/s', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_surf_PM25_sulf = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11 /), & + 'surf_PM25_sulf', 'Sulfate PM2.5 concentration at the surface', 'ug/m3', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_p_tropopause = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11 /), & + 'p_tropopause', 'Tropopause pressure', 'Pa', (/ ('', i=1, 9) /)) + TYPE(ctrl_out), SAVE :: o_sfluxaer = ctrl_out((/ 1, 10, 10, 10, 10, 10, 11, 11, 11 /), & + 'sflux', 'Ground sedimentation flux of strat. particles', 'kg(S)/m2/s', (/ ('', i=1, 9) /)) +#endif !!!!!!!!!!!!!!!!!!!!!! 3D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Index: LMDZ5/trunk/libf/phylmd/phys_output_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_output_mod.F90 (revision 2689) +++ LMDZ5/trunk/libf/phylmd/phys_output_mod.F90 (revision 2690) @@ -427,5 +427,5 @@ DO iq=nqo+1,nqtot iiq=niadv(iq) - o_trac(iq-nqo) = ctrl_out((/ 4, 5, 5, 5, 10, 10, 11, 11, 11 /), & + o_trac(iq-nqo) = ctrl_out((/ 1, 5, 5, 5, 10, 10, 11, 11, 11 /), & tname(iiq),'Tracer '//ttext(iiq), "-", & (/ '', '', '', '', '', '', '', '', '' /)) @@ -500,5 +500,5 @@ (/ '', '', '', '', '', '', '', '', '' /)) - o_trac_cum(iq-nqo) = ctrl_out((/ 3, 4, 10, 10, 10, 10, 11, 11, 11 /), & + o_trac_cum(iq-nqo) = ctrl_out((/ 1, 4, 10, 10, 10, 10, 11, 11, 11 /), & 'cum'//tname(iiq),& 'Cumulated tracer '//ttext(iiq), "-", & Index: LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 2689) +++ LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 2690) @@ -182,4 +182,12 @@ o_alt_tropo +#ifdef CPP_StratAer + USE phys_output_ctrlout_mod, only: & + o_sulf_convert, o_sulf_nucl, o_sulf_cond_evap, o_ocs_convert, & + o_sfluxaer, o_R2SO4, o_OCS_lifetime, o_SO2_lifetime, & + o_OCS_backgr_tend, o_SO2_backgr_tend, o_sulf_dep_dry, o_sulf_dep_wet, & + o_surf_PM25_sulf, o_ext_strat_550, o_tau_strat_550, & + o_p_tropopause, o_vsed_aer, o_tau_strat_1020, o_ext_strat_1020, o_f_r_wet +#endif USE phys_state_var_mod, only: pctsrf, paire_ter, rain_fall, snow_fall, & @@ -273,4 +281,13 @@ ep, epmax_diag ! epmax_cape +#ifdef CPP_StratAer + USE phys_local_var_mod, only: & + sulf_convert, sulf_nucl, sulf_cond_evap, ocs_convert, & + sfluxaer, R2SO4, OCS_lifetime, SO2_lifetime, & + OCS_backgr_tend, SO2_backgr_tend, sulf_dep_dry, sulf_dep_wet, & + surf_PM25_sulf, tau_strat_550, p_tropopause, tausum_strat, & + vsed_aer, tau_strat_1020, f_r_wet +#endif + USE phys_output_var_mod, only: vars_defined, snow_o, zfra_o, bils_diss, & bils_ec,bils_ech, bils_tke, bils_kinetic, bils_latent, bils_enthalp, & @@ -292,6 +309,4 @@ alt_tropo - - USE ocean_slab_mod, only: nslay, tslab, slab_bils, slab_bilg, tice, & seaice, slab_ekman,slab_hdiff, dt_ekman, dt_hdiff @@ -370,5 +385,5 @@ CALL set_itau_iophy(itau_w) - IF(.NOT.vars_defined) THEN + IF (.NOT.vars_defined) THEN iinitend = 2 ELSE @@ -381,10 +396,10 @@ !$OMP MASTER IF (vars_defined) THEN - if (prt_level >= 10) then + IF (prt_level >= 10) then write(lunout,*)"phys_output_write: call xios_update_calendar, itau_w=",itau_w - endif + ENDIF ! CALL xios_update_calendar(itau_w) CALL xios_update_calendar(itap) - END IF + ENDIF !$OMP END MASTER !$OMP BARRIER @@ -395,10 +410,10 @@ zx_tmp_fi2d = cell_area - if (is_north_pole_phy) then + IF (is_north_pole_phy) then zx_tmp_fi2d(1) = cell_area(1)/nbp_lon - endif - if (is_south_pole_phy) then + ENDIF + IF (is_south_pole_phy) then zx_tmp_fi2d(klon) = cell_area(klon)/nbp_lon - endif + ENDIf CALL histwrite_phy(o_aire, zx_tmp_fi2d) @@ -779,14 +794,14 @@ CALL histwrite_phy(o_uq, uq) CALL histwrite_phy(o_vq, vq) - IF(iflag_con.GE.3) THEN ! sb + IF (iflag_con.GE.3) THEN ! sb CALL histwrite_phy(o_cape, cape) CALL histwrite_phy(o_pbase, ema_pcb) CALL histwrite_phy(o_ptop, ema_pct) CALL histwrite_phy(o_fbase, ema_cbmf) - if (iflag_con /= 30) then + IF (iflag_con /= 30) THEN CALL histwrite_phy(o_plcl, plcl) CALL histwrite_phy(o_plfc, plfc) CALL histwrite_phy(o_wbeff, wbeff) - end if + ENDIF CALL histwrite_phy(o_cape_max, cape) @@ -799,5 +814,5 @@ CALL histwrite_phy(o_ftime_con, zx_tmp_fi2d) IF (vars_defined) THEN - IF(iflag_thermals>=1)THEN + IF (iflag_thermals>=1)THEN zx_tmp_fi3d=dnwd+dnwd0+upwd+fm_therm(:,1:klev) ELSE @@ -850,5 +865,5 @@ DO k=1, nlevSTD bb2=clevSTD(k) - IF(bb2.EQ."850".OR.bb2.EQ."700".OR. & + IF (bb2.EQ."850".OR.bb2.EQ."700".OR. & bb2.EQ."500".OR.bb2.EQ."200".OR. & bb2.EQ."100".OR. & @@ -871,5 +886,5 @@ #endif #ifdef CPP_XIOS - IF(ok_all_xml) THEN + IF (ok_all_xml) THEN !XIOS CALL xios_get_field_attr("u850",default_value=missing_val) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -877,5 +892,5 @@ DO k=1, nlevSTD bb2=clevSTD(k) - IF(bb2.EQ."850".OR.bb2.EQ."700".OR. & + IF (bb2.EQ."850".OR.bb2.EQ."700".OR. & bb2.EQ."500".OR.bb2.EQ."200".OR. & bb2.EQ."100".OR. & @@ -994,10 +1009,10 @@ ELSE CALL histwrite_phy(o_tslab, tslab(:,1:nslay)) - END IF + ENDIF IF (version_ocean=='sicINT') THEN CALL histwrite_phy(o_slab_bilg, slab_bilg) CALL histwrite_phy(o_slab_tice, tice) CALL histwrite_phy(o_slab_sic, seaice) - END IF + ENDIF IF (slab_hdiff) THEN IF (nslay.EQ.1) THEN @@ -1006,6 +1021,6 @@ ELSE CALL histwrite_phy(o_slab_hdiff, dt_hdiff(:,1:nslay)) - END IF - END IF + ENDIF + ENDIF IF (slab_ekman.GT.0) THEN IF (nslay.EQ.1) THEN @@ -1014,6 +1029,6 @@ ELSE CALL histwrite_phy(o_slab_ekman, dt_ekman(:,1:nslay)) - END IF - END IF + ENDIF + ENDIF ENDIF !type_ocean == force/slab CALL histwrite_phy(o_weakinv, weak_inversion) @@ -1094,16 +1109,36 @@ ENDIF IF (flag_aerosol.GT.0.OR.flag_aerosol_strat.GT.0) THEN -! DO naero = 1, naero_spc -!--correction mini bug OB DO naero = 1, naero_tot - CALL histwrite_phy(o_tausumaero(naero), & - tausum_aero(:,2,naero) ) + CALL histwrite_phy(o_tausumaero(naero),tausum_aero(:,2,naero)) END DO ENDIF IF (flag_aerosol_strat.GT.0) THEN - CALL histwrite_phy(o_tausumaero_lw, & - tausum_aero(:,6,id_STRAT_phy) ) - ENDIF - ENDIF + CALL histwrite_phy(o_tausumaero_lw,tausum_aero(:,6,id_STRAT_phy)) + ENDIF + ENDIF +#ifdef CPP_StratAer + IF (type_trac=='coag') THEN + CALL histwrite_phy(o_sulf_convert, sulf_convert) + CALL histwrite_phy(o_sulf_nucl, sulf_nucl) + CALL histwrite_phy(o_sulf_cond_evap, sulf_cond_evap) + CALL histwrite_phy(o_ocs_convert, ocs_convert) + CALL histwrite_phy(o_R2SO4, R2SO4) + CALL histwrite_phy(o_OCS_lifetime, OCS_lifetime) + CALL histwrite_phy(o_SO2_lifetime, SO2_lifetime) + CALL histwrite_phy(o_OCS_backgr_tend, OCS_backgr_tend) + CALL histwrite_phy(o_SO2_backgr_tend, SO2_backgr_tend) + CALL histwrite_phy(o_sulf_dep_dry, sulf_dep_dry) + CALL histwrite_phy(o_sulf_dep_wet, sulf_dep_wet) + CALL histwrite_phy(o_surf_PM25_sulf, surf_PM25_sulf) + CALL histwrite_phy(o_p_tropopause, p_tropopause) + CALL histwrite_phy(o_sfluxaer, sfluxaer) + CALL histwrite_phy(o_vsed_aer, vsed_aer) + CALL histwrite_phy(o_f_r_wet, f_r_wet) + CALL histwrite_phy(o_ext_strat_550, tau_strat_550) + CALL histwrite_phy(o_ext_strat_1020, tau_strat_1020) + CALL histwrite_phy(o_tau_strat_550, tausum_strat(:,1)) + CALL histwrite_phy(o_tau_strat_1020, tausum_strat(:,2)) + ENDIF +#endif IF (ok_ade) THEN CALL histwrite_phy(o_topswad, topswad_aero*swradcorr) @@ -1116,5 +1151,5 @@ CALL histwrite_phy(o_sollwad0, sollwad0_aero) !====MS forcing diagnostics - if (new_aod) then + IF (new_aod) THEN zx_tmp_fi2d(:)=topsw_aero(:,1)*swradcorr(:) CALL histwrite_phy(o_swtoaas_nat,zx_tmp_fi2d) @@ -1135,5 +1170,5 @@ CALL histwrite_phy(o_swsrfcs_ant,zx_tmp_fi2d) !cf - if (.not. aerosol_couple) then + IF (.not. aerosol_couple) THEN zx_tmp_fi2d(:)=topswcf_aero(:,1)*swradcorr(:) CALL histwrite_phy(o_swtoacf_nat,zx_tmp_fi2d) @@ -1148,6 +1183,6 @@ zx_tmp_fi2d(:)=solswcf_aero(:,3)*swradcorr(:) CALL histwrite_phy(o_swsrfcf_zero,zx_tmp_fi2d) - endif - endif ! new_aod + ENDIF + ENDIF ! new_aod !====MS forcing diagnostics ENDIF @@ -1261,13 +1296,13 @@ CALL histwrite_phy(o_alb2, albsol2) !FH Sorties pour la couche limite - if (iflag_pbl>1) then + IF (iflag_pbl>1) THEN zx_tmp_fi3d=0. IF (vars_defined) THEN - do nsrf=1,nbsrf - do k=1,klev + DO nsrf=1,nbsrf + DO k=1,klev zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) & +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf) - enddo - enddo + ENDDO + ENDDO ENDIF CALL histwrite_phy(o_tke, zx_tmp_fi3d) @@ -1325,8 +1360,8 @@ CALL histwrite_phy(o_dqcon2d, zx_tmp_fi2d) - IF(iflag_thermals.EQ.0) THEN + IF (iflag_thermals.EQ.0) THEN IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_con(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_tnhusc, zx_tmp_fi3d) - ELSE IF(iflag_thermals.GE.1.AND.iflag_wake.EQ.1) THEN + ELSE IF (iflag_thermals.GE.1.AND.iflag_wake.EQ.1) THEN IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)=d_q_con(1:klon,1:klev)/pdtphys + & @@ -1350,14 +1385,14 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Sorties specifiques a la separation thermiques/non thermiques - if (iflag_thermals>=1) then - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lscth(1:klon,1:klev)/pdtphys + IF (iflag_thermals>=1) THEN + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lscth(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtlscth, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lscst(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lscst(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtlscst, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_lscth(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_lscth(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dqlscth, zx_tmp_fi3d) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d) CALL histwrite_phy(o_dqlscth2d, zx_tmp_fi2d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_lscst(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_lscst(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dqlscst, zx_tmp_fi3d) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d) @@ -1366,13 +1401,13 @@ CALL histwrite_phy(o_plulst, plul_st) IF (vars_defined) THEN - do k=1,klev - do i=1,klon - if (ptconvth(i,k)) then + DO k=1,klev + DO i=1,klon + IF (ptconvth(i,k)) THEN zx_tmp_fi3d(i,k)=1. - else + ELSE zx_tmp_fi3d(i,k)=0. - endif - enddo - enddo + ENDIF + ENDDO + ENDDO ENDIF CALL histwrite_phy(o_ptconvth, zx_tmp_fi3d) @@ -1383,9 +1418,9 @@ ENDIF CALL histwrite_phy(o_lmaxth, zx_tmp_fi2d) - endif ! iflag_thermals>=1 + ENDIF ! iflag_thermals>=1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_vdf(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtvdf, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_diss(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_diss(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtdis, zx_tmp_fi3d) IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_vdf(1:klon,1:klev)/pdtphys @@ -1437,40 +1472,40 @@ CALL histwrite_phy(o_dqthe2d, zx_tmp_fi2d) ENDIF !iflag_thermals - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtajs, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dqajs, zx_tmp_fi3d) CALL water_int(klon,klev,zx_tmp_fi3d,zmasse,zx_tmp_fi2d) CALL histwrite_phy(o_dqajs2d, zx_tmp_fi2d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_swr(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_swr(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtswr, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_sw0(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_sw0(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtsw0, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lwr(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lwr(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtlwr, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lw0(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lw0(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtlw0, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_ec(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_ec(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtec, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_u_vdf(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_u_vdf(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_duvdf, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_v_vdf(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_v_vdf(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dvvdf, zx_tmp_fi3d) IF (ok_orodr) THEN - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_u_oro(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_u_oro(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_duoro, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_v_oro(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_v_oro(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dvoro, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_oro(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_oro(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtoro, zx_tmp_fi3d) ENDIF IF (ok_orolf) THEN - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_u_lif(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_u_lif(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dulif, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_v_lif(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_v_lif(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dvlif, zx_tmp_fi3d) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lif(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lif(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dtlif, zx_tmp_fi3d) ENDIF @@ -1479,11 +1514,11 @@ CALL histwrite_phy(o_du_gwd_hines, du_gwd_hines/pdtphys) CALL histwrite_phy(o_dv_gwd_hines, dv_gwd_hines/pdtphys) - IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_hin(1:klon,1:klev)/pdtphys + IF (vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_hin(1:klon,1:klev)/pdtphys CALL histwrite_phy(o_dthin, zx_tmp_fi3d) CALL histwrite_phy(o_ustr_gwd_hines, zustr_gwd_hines) CALL histwrite_phy(o_vstr_gwd_hines, zvstr_gwd_hines) - end IF - - if (.not. ok_hines .and. ok_gwd_rando) then + ENDIF + + IF (.not. ok_hines .and. ok_gwd_rando) THEN CALL histwrite_phy(o_du_gwd_front, du_gwd_front / pdtphys) CALL histwrite_phy(o_dv_gwd_front, dv_gwd_front / pdtphys) @@ -1492,5 +1527,5 @@ ENDIF - IF (ok_gwd_rando) then + IF (ok_gwd_rando) THEN CALL histwrite_phy(o_du_gwd_rando, du_gwd_rando / pdtphys) CALL histwrite_phy(o_dv_gwd_rando, dv_gwd_rando / pdtphys) @@ -1499,7 +1534,7 @@ CALL histwrite_phy(o_east_gwstress, east_gwstress ) CALL histwrite_phy(o_west_gwstress, west_gwstress ) - end IF - - IF (ok_qch4) then + ENDIF + + IF (ok_qch4) THEN CALL histwrite_phy(o_dqch4, d_q_ch4 / pdtphys) ENDIF @@ -1527,5 +1562,5 @@ CALL histwrite_phy(o_rldcs, lwdn0) - IF(vars_defined) THEN + IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)=d_t(1:klon,1:klev)+ & d_t_dyn(1:klon,1:klev) @@ -1533,10 +1568,10 @@ CALL histwrite_phy(o_tnt, zx_tmp_fi3d) - IF(vars_defined) THEN + IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)=d_t_swr(1:klon,1:klev)/pdtphys + & d_t_lwr(1:klon,1:klev)/pdtphys ENDIF CALL histwrite_phy(o_tntr, zx_tmp_fi3d) - IF(vars_defined) THEN + IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)= (d_t_lsc(1:klon,1:klev)+ & d_t_eva(1:klon,1:klev)+ & @@ -1544,10 +1579,10 @@ ENDIF CALL histwrite_phy(o_tntscpbl, zx_tmp_fi3d) - IF(vars_defined) THEN + IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)=d_qx(1:klon,1:klev,ivap)+ & d_q_dyn(1:klon,1:klev) ENDIF CALL histwrite_phy(o_tnhus, zx_tmp_fi3d) - IF(vars_defined) THEN + IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)=d_q_lsc(1:klon,1:klev)/pdtphys+ & d_q_eva(1:klon,1:klev)/pdtphys @@ -1555,36 +1590,36 @@ CALL histwrite_phy(o_tnhusscpbl, zx_tmp_fi3d) CALL histwrite_phy(o_evu, coefm(:,:,is_ave)) - IF(vars_defined) THEN + IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)=q_seri(1:klon,1:klev)+ & ql_seri(1:klon,1:klev) ENDIF CALL histwrite_phy(o_h2o, zx_tmp_fi3d) - if (iflag_con >= 3) then - IF(vars_defined) THEN + IF (iflag_con >= 3) THEN + IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)=-1 * (dnwd(1:klon,1:klev)+ & dnwd0(1:klon,1:klev)) ENDIF CALL histwrite_phy(o_mcd, zx_tmp_fi3d) - IF(vars_defined) THEN + IF (vars_defined) THEN zx_tmp_fi3d(1:klon,1:klev)=upwd(1:klon,1:klev) + & dnwd(1:klon,1:klev)+ dnwd0(1:klon,1:klev) ENDIF CALL histwrite_phy(o_dmc, zx_tmp_fi3d) - else if (iflag_con == 2) then + ELSE IF (iflag_con == 2) THEN CALL histwrite_phy(o_mcd, pmfd) CALL histwrite_phy(o_dmc, pmfu + pmfd) - end if + ENDIF CALL histwrite_phy(o_ref_liq, ref_liq) CALL histwrite_phy(o_ref_ice, ref_ice) - if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. & + IF (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. & RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. & RCFC12_per.NE.RCFC12_act) THEN - IF(vars_defined) zx_tmp_fi2d(:) = swupp(:,klevp1)*swradcorr(:) + IF (vars_defined) zx_tmp_fi2d(:) = swupp(:,klevp1)*swradcorr(:) CALL histwrite_phy(o_rsut4co2, zx_tmp_fi2d) - IF(vars_defined) zx_tmp_fi2d(:) = lwupp(:,klevp1) + IF (vars_defined) zx_tmp_fi2d(:) = lwupp(:,klevp1) CALL histwrite_phy(o_rlut4co2, zx_tmp_fi2d) - IF(vars_defined) zx_tmp_fi2d(:) = swup0p(:,klevp1)*swradcorr(:) + IF (vars_defined) zx_tmp_fi2d(:) = swup0p(:,klevp1)*swradcorr(:) CALL histwrite_phy(o_rsutcs4co2, zx_tmp_fi2d) - IF(vars_defined) zx_tmp_fi2d(:) = lwup0p(:,klevp1) + IF (vars_defined) zx_tmp_fi2d(:) = lwup0p(:,klevp1) CALL histwrite_phy(o_rlutcs4co2, zx_tmp_fi2d) DO k=1, klevp1 @@ -1626,9 +1661,9 @@ CALL histwrite_phy(o_va,vwriteSTD(:,:,iff-6),iff) CALL histwrite_phy(o_wap,wwriteSTD(:,:,iff-6),iff) - IF(vars_defined) THEN + IF (vars_defined) THEN DO k=1, nlevSTD DO i=1, klon - IF(tnondef(i,k,iff-6).NE.missing_val) THEN - IF(freq_outNMC(iff-6).LT.0) THEN + IF (tnondef(i,k,iff-6).NE.missing_val) THEN + IF (freq_outNMC(iff-6).LT.0) THEN freq_moyNMC(iff-6)=(mth_len*un_jour)/freq_calNMC(iff-6) ELSE @@ -1643,8 +1678,8 @@ ENDIF CALL histwrite_phy(o_psbg,zx_tmp_fi3d_STD,iff) - IF(vars_defined) THEN + IF (vars_defined) THEN DO k=1, nlevSTD DO i=1, klon - IF(O3sumSTD(i,k,iff-6).NE.missing_val) THEN + IF (O3sumSTD(i,k,iff-6).NE.missing_val) THEN zx_tmp_fi3d_STD(i,k) = O3sumSTD(i,k,iff-6) * 1.e+9 ELSE @@ -1655,9 +1690,9 @@ ENDIF CALL histwrite_phy(o_tro3,zx_tmp_fi3d_STD,iff) - if (read_climoz == 2) THEN - IF(vars_defined) THEN + IF (read_climoz == 2) THEN + IF (vars_defined) THEN DO k=1, nlevSTD DO i=1, klon - IF(O3daysumSTD(i,k,iff-6).NE.missing_val) THEN + IF (O3daysumSTD(i,k,iff-6).NE.missing_val) THEN zx_tmp_fi3d_STD(i,k) = O3daysumSTD(i,k,iff-6) * 1.e+9 ELSE @@ -1683,5 +1718,5 @@ #endif #ifdef CPP_XIOS - IF(ok_all_xml) THEN + IF (ok_all_xml) THEN ! DO iff=7, nfiles @@ -1694,9 +1729,9 @@ CALL histwrite_phy(o_va,vlevSTD(:,:)) CALL histwrite_phy(o_wap,wlevSTD(:,:)) -! IF(vars_defined) THEN +! IF (vars_defined) THEN ! DO k=1, nlevSTD ! DO i=1, klon -! IF(tnondef(i,k,3).NE.missing_val) THEN -! IF(freq_outNMC(iff-6).LT.0) THEN +! IF (tnondef(i,k,3).NE.missing_val) THEN +! IF (freq_outNMC(iff-6).LT.0) THEN ! freq_moyNMC(iff-6)=(mth_len*un_jour)/freq_calNMC(iff-6) ! ELSE @@ -1711,8 +1746,8 @@ ! ENDIF ! CALL histwrite_phy(o_psbg,zx_tmp_fi3d_STD) - IF(vars_defined) THEN + IF (vars_defined) THEN DO k=1, nlevSTD DO i=1, klon - IF(O3STD(i,k).NE.missing_val) THEN + IF (O3STD(i,k).NE.missing_val) THEN zx_tmp_fi3d_STD(i,k) = O3STD(i,k) * 1.e+9 ELSE @@ -1723,9 +1758,9 @@ ENDIF CALL histwrite_phy(o_tro3,zx_tmp_fi3d_STD) - if (read_climoz == 2) THEN - IF(vars_defined) THEN + IF (read_climoz == 2) THEN + IF (vars_defined) THEN DO k=1, nlevSTD DO i=1, klon - IF(O3daySTD(i,k).NE.missing_val) THEN + IF (O3daySTD(i,k).NE.missing_val) THEN zx_tmp_fi3d_STD(i,k) = O3daySTD(i,k) * 1.e+9 ELSE @@ -1736,5 +1771,5 @@ ENDIF CALL histwrite_phy(o_tro3_daylight,zx_tmp_fi3d_STD) - endif + ENDIF CALL histwrite_phy(o_uxv,uvSTD(:,:)) CALL histwrite_phy(o_vxq,vqSTD(:,:)) @@ -1751,11 +1786,7 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (nqtot.GE.nqo+1) THEN - DO iq=nqo+1,nqtot - IF (type_trac == 'lmdz' .OR. type_trac == 'repr') THEN - -!jyg< -!! CALL histwrite_phy(o_trac(iq-nqo), qx(:,:,iq)) + DO iq=nqo+1,nqtot + IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag') THEN CALL histwrite_phy(o_trac(iq-nqo), tr_seri(:,:,iq-nqo)) -!>jyg CALL histwrite_phy(o_dtr_vdf(iq-nqo),d_tr_cl(:,:,iq-nqo)) CALL histwrite_phy(o_dtr_the(iq-nqo),d_tr_th(:,:,iq-nqo)) @@ -1772,18 +1803,15 @@ CALL histwrite_phy(o_dtr_uscav(iq-nqo),d_tr_uscav(:,:,iq-nqo)) zx_tmp_fi2d=0. - IF(vars_defined) THEN + IF (vars_defined) THEN DO k=1,klev -!jyg< -!! zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*qx(:,k,iq) zx_tmp_fi2d(:)=zx_tmp_fi2d(:)+zmasse(:,k)*tr_seri(:,k,iq-nqo) -!>jyg ENDDO ENDIF CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d) - endif - ENDDO - ENDIF - - IF(.NOT.vars_defined) THEN + ENDIF + ENDDO + ENDIF + + IF (.NOT.vars_defined) THEN !$OMP MASTER #ifndef CPP_IOIPSL_NO_OUTPUT @@ -1801,14 +1829,12 @@ CALL wxios_closedef() #endif - !$OMP END MASTER !$OMP BARRIER vars_defined = .TRUE. - - END IF - - END DO - - IF(vars_defined) THEN + ENDIF + + ENDDO + + IF (vars_defined) THEN ! On synchronise les fichiers pour IOIPSL #ifndef CPP_IOIPSL_NO_OUTPUT @@ -1823,5 +1849,4 @@ ENDIF - END SUBROUTINE phys_output_write Index: LMDZ5/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/physiq_mod.F90 (revision 2689) +++ LMDZ5/trunk/libf/phylmd/physiq_mod.F90 (revision 2690) @@ -745,7 +745,5 @@ REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s) REAL conv_t(klon,klev) ! convergence de la temperature(K/s) - - ! - ! REAL zxsnow(klon) + ! REAL zxsnow_dummy(klon) REAL zsav_tsol(klon) @@ -763,5 +761,5 @@ real zqsat(klon,klev) ! - INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq + INTEGER i, k, iq, j, nsrf, ll, l ! REAL t_coup @@ -884,5 +882,4 @@ !IM 141004 END !IM 190504 BEG - INTEGER ij ! INTEGER imp1jmp1 ! PARAMETER(imp1jmp1=(iim+1)*jjmp1) @@ -893,6 +890,4 @@ LOGICAL ok_msk REAL msk(klon) - !IM - REAL airetot, pi !ym A voir plus tard !ym REAL zm_wo(jjmp1, klev) @@ -932,5 +927,4 @@ !$OMP THREADPRIVATE(ok_sync) real date0 - integer idayref ! essai writephys @@ -953,7 +947,4 @@ DATA ip_ebil/0/ !$OMP THREADPRIVATE(ip_ebil) - INTEGER if_ebil ! level for energy conserv. dignostics - SAVE if_ebil - !$OMP THREADPRIVATE(if_ebil) REAL q2m(klon,nbsrf) ! humidite a 2m @@ -3483,4 +3474,7 @@ ELSE #ifdef CPP_RRTM +#ifndef CPP_StratAer + !--prescribed strat aerosols + !--only in the case of non-interactive strat aerosols IF (flag_aerosol_strat.EQ.1) THEN CALL readaerosolstrato1_rrtm(debut) @@ -3492,4 +3486,5 @@ CALL abort_physic(modname,abort_message,1) ENDIF +#endif #else abort_message='You should compile with -rrtm if running ' & @@ -3499,4 +3494,11 @@ ENDIF ENDIF +! +#ifdef CPP_RRTM +#ifdef CPP_StratAer + !--interactive strat aerosols + CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut) +#endif +#endif !--fin STRAT AEROSOL ! Index: LMDZ5/trunk/libf/phylmd/phytrac_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phytrac_mod.F90 (revision 2689) +++ LMDZ5/trunk/libf/phylmd/phytrac_mod.F90 (revision 2690) @@ -97,8 +97,14 @@ USE tracreprobus_mod USE indice_sol_mod - USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root USE print_control_mod, ONLY: lunout USE aero_mod, ONLY : naero_grp + +#ifdef CPP_StratAer + USE traccoag_mod + USE phys_local_var_mod, ONLY: mdw, sulf_dep_dry, sulf_dep_wet + USE infotrac, ONLY: nbtr_sulgas, id_SO2_strat, id_H2SO4_strat + USE aerophys +#endif IMPLICIT NONE @@ -208,5 +214,4 @@ !-------------- ! - ! REAL,DIMENSION(:),INTENT(IN) :: cdragh ! (klon) coeff drag pour T et Q REAL,DIMENSION(:,:),INTENT(IN) :: coefh ! (klon,klev) coeff melange CL (m**2/s) @@ -215,5 +220,4 @@ REAL,DIMENSION(:),INTENT(IN) :: yu1 ! (klon) vents au premier niveau REAL,DIMENSION(:),INTENT(IN) :: yv1 ! (klon) vents au premier niveau - ! !Lessivage: @@ -238,9 +242,12 @@ REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol) - +#ifdef CPP_StratAer + REAL,DIMENSION(klon) :: v_dep_dry !dry deposition velocity of stratospheric sulfate in m/s +#endif ! Output argument !---------------- REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] REAL,DIMENSION(klon,klev) :: sourceBE + !======================================================================================= ! -- LOCAL VARIABLES -- @@ -267,5 +274,4 @@ INTEGER :: itau_w ! pas de temps ecriture = nstep + itau_phy LOGICAL,PARAMETER :: ok_sync=.TRUE. - ! ! Nature du traceur @@ -456,4 +462,15 @@ CASE('repr') source(:,:)=0. +#ifdef CPP_StratAer + CASE('coag') + source(:,:)=0. + DO it= 1, nbtr_sulgas + aerosol(it)=.FALSE. + IF (it==id_H2SO4_strat) aerosol(it)=.TRUE. + ENDDO + DO it= nbtr_sulgas+1, nbtr + aerosol(it)=.TRUE. + ENDDO +#endif END SELECT @@ -504,4 +521,17 @@ !--for now we do not scavenge in cvltr flag_cvltr(it)=.false. + +#ifdef CPP_StratAer + CASE('coag') + IF (convscav.and.aerosol(it)) THEN + flag_cvltr(it)=.true. + ccntrAA(it) =ccntrAA_in + ccntrENV(it)=ccntrENV_in + coefcoli(it)=coefcoli_in + ELSE + flag_cvltr(it)=.false. + ENDIF +#endif + END SELECT ENDDO @@ -572,7 +602,7 @@ ! Appel fait en fin de phytrac pour avoir les emissions modifiees par ! la couche limite et la convection avant le calcul de la chimie + CASE('repr') ! -- CHIMIE REPROBUS -- - CALL tracreprobus(pdtphys, gmtime, debutphy, julien, & presnivs, xlat, xlon, pphis, pphi, & @@ -580,4 +610,13 @@ tr_seri) +#ifdef CPP_StratAer + CASE('coag') + ! --STRATOSPHERIC AER IN THE STRAT -- + CALL traccoag(pdtphys, gmtime, debutphy, julien, & + presnivs, xlat, xlon, pphis, pphi, & + t_seri, pplay, paprs, sh, rh , & + tr_seri) +#endif + END SELECT !====================================================================== @@ -591,6 +630,6 @@ IF (iflag_con.LT.2) THEN !--pas de transport convectif - d_tr_cv(:,:,it)=0. + ELSE IF (iflag_con.EQ.2) THEN !--ancien transport convectif de Tiedtke @@ -648,4 +687,28 @@ END DO ! nbtr + +#ifdef CPP_StratAer + IF (type_trac=='coag') THEN + ! initialize wet deposition flux of sulfur + sulf_dep_wet(:)=0.0 + ! compute wet deposition flux of sulfur (sum over gases and particles) + ! and convert to kg(S)/m2/s + DO i = 1, klon + DO k = 1, klev + DO it = 1, nbtr + !do not include SO2 because most of it comes trom the troposphere + IF (it==id_H2SO4_strat) THEN + sulf_dep_wet(i)=sulf_dep_wet(i)+d_tr_cv(i,k,it)*(mSatom/mH2SO4mol) & + & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys + ELSEIF (it.GT.nbtr_sulgas) THEN + sulf_dep_wet(i)=sulf_dep_wet(i)+d_tr_cv(i,k,it)*(mSatom/mH2SO4mol) & + & *dens_aer_dry*4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3 & + & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys + ENDIF + ENDDO + ENDDO + ENDDO + ENDIF +#endif END IF ! convection @@ -692,4 +755,27 @@ ! Injection during BL mixing ! +#ifdef CPP_StratAer + IF (type_trac=='coag') THEN + + ! initialize dry deposition flux of sulfur + sulf_dep_dry(:)=0.0 + + ! compute dry deposition velocity as function of surface type (numbers + ! from IPSL note 23, 2002) + v_dep_dry(:) = pctsrf(:,is_ter) * 2.5e-3 & + & + pctsrf(:,is_oce) * 0.5e-3 & + & + pctsrf(:,is_lic) * 2.5e-3 & + & + pctsrf(:,is_sic) * 2.5e-3 + + ! compute surface dry deposition flux + zrho(:,1)=pplay(:,1)/t_seri(:,1)/RD + + DO it=1, nbtr + source(:,it) = - v_dep_dry(:) * tr_seri(:,1,it) * zrho(:,1) + ENDDO + + ENDIF +#endif + DO it=1, nbtr ! @@ -703,7 +789,19 @@ tr_seri(:,:,it)=tr_seri(:,:,it)+d_tr_cl(:,:,it) ! - END IF +#ifdef CPP_StratAer + IF (type_trac=='coag') THEN + ! compute dry deposition flux of sulfur (sum over gases and particles) + IF (it==id_H2SO4_strat) THEN + sulf_dep_dry(:)=sulf_dep_dry(:)-source(:,it)*(mSatom/mH2SO4mol) + ELSEIF (it.GT.nbtr_sulgas) THEN + sulf_dep_dry(:)=sulf_dep_dry(:)-source(:,it)*(mSatom/mH2SO4mol)*dens_aer_dry & + & *4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3 + ENDIF + ENDIF +#endif + ! + ENDIF ! - END DO + ENDDO ! ELSE IF (iflag_vdf_trac==0) THEN @@ -720,5 +818,4 @@ ! ! Nothing happens - ! d_tr_cl=0. ! @@ -772,4 +869,26 @@ END DO !tr + +#ifdef CPP_StratAer + IF (type_trac=='coag') THEN + ! compute wet deposition flux of sulfur (sum over gases and + ! particles) and convert to kg(S)/m2/s + ! adding contribution of d_tr_ls to d_tr_cv (above) + DO i = 1, klon + DO k = 1, klev + DO it = 1, nbtr + IF (it==id_H2SO4_strat) THEN + sulf_dep_wet(i)=sulf_dep_wet(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol) & + & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys + ELSEIF (it.GT.nbtr_sulgas) THEN + sulf_dep_wet(i)=sulf_dep_wet(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol) & + & *dens_aer_dry*4./3.*RPI*(mdw(it-nbtr_sulgas)/2.)**3 & + & *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys + ENDIF + ENDDO + ENDDO + ENDDO + ENDIF +#endif ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl Index: LMDZ5/trunk/makelmdz =================================================================== --- LMDZ5/trunk/makelmdz (revision 2689) +++ LMDZ5/trunk/makelmdz (revision 2690) @@ -29,4 +29,5 @@ rrtm=false dust=false +strataer=false full="" @@ -113,4 +114,5 @@ [-rrtm true/false] : compile with/without rrtm package (default: false) [-dust true/false] : compile with/without the dust package from Boucher et al. (default: false) +[-strataer true/false] : compile with/without the strat aer package from Boucher et al. (default: false) [-parallel none/mpi/omp/mpi_omp] : parallelism (default: none) : mpi, openmp or mixted mpi_openmp [-g GRI] : grid configuration in dyn3d/GRI_xy.h (default: reg, inclues a zoom) @@ -180,4 +182,7 @@ dust="$2" ; shift ; shift ;; + "-strataer") + strataer="$2" ; shift ; shift ;; + "-mem") paramem="mem" ; shift ;; @@ -459,4 +464,10 @@ CPP_KEY="$CPP_KEY CPP_Dust" src_dirs="$src_dirs phy${physique}/Dust" +fi + +if [[ "$strataer" == "true" ]] +then + CPP_KEY="$CPP_KEY CPP_StratAer" + src_dirs="$src_dirs phy${physique}/StratAer" fi Index: LMDZ5/trunk/makelmdz_fcm =================================================================== --- LMDZ5/trunk/makelmdz_fcm (revision 2689) +++ LMDZ5/trunk/makelmdz_fcm (revision 2690) @@ -26,4 +26,5 @@ rrtm=false dust=false +strataer=false chimie=false parallel=none @@ -48,4 +49,5 @@ RRTM_PATH=$LMDGCM/.void_dir DUST_PATH=$LMDGCM/.void_dir +STRATAER_PATH=$LMDGCM/.void_dir SISVAT_PATH=$LMDGCM/.void_dir COSP_PATH=$LMDGCM/.void_dir @@ -90,4 +92,5 @@ [-rrtm true/false] : compile with/without rrtm package (default: false) [-dust true/false] : compile with/without the dust package by Boucher and co (default: false) +[-strataer true/false] : compile with/without the strat aer package by Boucher and co (default: false) [-parallel none/mpi/omp/mpi_omp] : parallelism (default: none) : mpi, openmp or mixted mpi_openmp [-g GRI] : grid configuration in dyn3d/GRI_xy.h (default: reg, inclues a zoom) @@ -144,4 +147,7 @@ "-dust") dust="$2" ; shift ; shift ;; + + "-strataer") + strataer="$2" ; shift ; shift ;; "-chimie") @@ -363,4 +369,10 @@ fi +if [[ "$strataer" == "true" ]] +then + CPP_KEY="$CPP_KEY CPP_StratAer" + STRATAER_PATH="$LIBFGCM/%PHYS/StratAer" +fi + if [[ $io == ioipsl ]] then @@ -593,4 +605,5 @@ echo "%RRTM $RRTM_PATH" >> $config_fcm echo "%DUST $DUST_PATH" >> $config_fcm +echo "%STRATAER $STRATAER_PATH" >> $config_fcm echo "%SISVAT $SISVAT_PATH" >> $config_fcm echo "%COSP $COSP_PATH" >> $config_fcm