subroutine hazecloud(ngrid,nlayer,nq,ptimestep,zday, & pplay,pplev,zzlay,pq,pdq,pdist_sol,mu0,pfluxuv,zdqhaze, & zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin) !zdqhaze_col) use comgeomfi_h use comcstfi_mod, only: pi, g use tracer_h, only: igcm_haze, igcm_ch4_gas, igcm_n2, igcm_prec_haze, noms, mmol use geometry_mod, only: longitude, latitude ! in radians use callkeys_mod, only: hazeconservch4, haze_ch4proffix, diurnal, tcon_ch4, k_ch4, & ncratio_ch4,triton,global1d use datafile_mod, only: datadir use write_output_mod, only: write_output implicit none !================================================================== ! Purpose ! ------- ! Production of haze in the atmosphere ! ! Inputs ! ------ ! ngrid Number of vertical columns ! nlayer Number of layers ! nq Number of tracers ! pplay(ngrid,nlayer) Pressure layers ! pplev(ngrid,nlayer+1) Pressure levels ! zzlay(ngrid,nlayer+1) mid layer altitude ! ptimestep Physical timestep ! zday, mu0 day clock, cos(incident flux angle) ! pq, pdq tracer MMR and tendency ! pdist_sol, declin Sun-planet distance, subsolar latitude ! pfluxuv Lyman alpha flux at Earth ! Outputs ! ------- ! zdqhaze Tendency on tracers MMR haze ! zdqphot_prec Photolysis rate for haze precursors ! zdqphot_ch4 Photolysis rate for CH4 ! zdqconv_prec Tendency for haze formation and precursor destruction ! Both ! ---- ! ! Authors ! ------- ! Tanguy Bertrand and Francois Forget (2014) ! !================================================================== !----------------------------------------------------------------------- ! Arguments INTEGER ngrid, nlayer, nq LOGICAL firstcall SAVE firstcall DATA firstcall/.true./ REAL ptimestep REAL zday REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1) REAL,INTENT(IN) :: zzlay(ngrid,nlayer) ! Mid-layer altitude REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) REAL,INTENT(IN) :: pdist_sol ! distance SUN-pluto in AU REAL,INTENT(IN) :: pfluxuv ! Lyman alpha flux at Earth at specific Ls (ph/cm2/s) REAL,INTENT(IN) :: mu0(ngrid) ! cosinus of solar incident flux REAL,INTENT(IN) :: declin ! Subsolar latitude REAL,INTENT(OUT) :: zdqhaze(ngrid,nlayer,nq) ! Final tendancy REAL,INTENT(OUT) :: zdqphot_ch4(ngrid,nlayer) ! tendancy due to photolysis of ch4 REAL,INTENT(OUT) :: zdqphot_prec(ngrid,nlayer) ! tendancy due to photolysis of ch4 REAL,INTENT(OUT) :: zdqconv_prec(ngrid,nlayer) ! tendancy due to conversion of ! prec_haze into haze !----------------------------------------------------------------------- ! Local variables INTEGER l,ig,iq REAL zq_ch4(ngrid,nlayer) REAL zq_prec(ngrid,nlayer) REAL tauch4(nlayer) REAL sigch4 ! aborption cross section ch4 cm-2 mol-1 REAL flym_sol_earth ! initial flux Earth ph.m-2.s-1 REAL flym_sol_pluto ! initial Lyman alpha SOLAR flux Pluto ph.m-2.s-1 REAL flym_ipm(ngrid) ! top of atm. Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 REAL fluxlym_sol(nlayer+1) ! Local Solar flux Lyman alpha ph.m-2.s-1 REAL fluxlym_ipm(nlayer+1) ! Local IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 REAL mu_ipm(ngrid) ! local Mean incident flux for IPM Lyman alpha photons REAL mu_sol(ngrid) ! local Mean incident flux for solar lyman alpha photons REAL gradflux(nlayer) ! gradient flux Lyman alpha ph.m-2.s-1 REAL kch4 ! fraction of Carbon from ch4 directly dissociated ! by prec_haze CHARACTER(len=10) :: tname REAL ncratio ! ration Nitrogen/Carbon observed in tholins REAL tcon ! Time constant: conversion in aerosol REAL avogadro ! avogadro constant integer Nfine,ifine parameter(Nfine=600) character(len=100) :: file_path real,save :: levdat(Nfine),vmrdat(Nfine) real :: vmrch4(ngrid,nlayer) ! vmr ch4 from vmrch4_proffix REAL longit REAL latit REAL valmin REAL valmax REAL valmin_dl REAL puis REAL dist REAL long2 ! Diagnostics (temporary) REAL fluxlym_sol_tot(ngrid) REAL fluxlym_ipm_tot(ngrid) !----------------------------------------------------------------------- ! Input parameters avogadro = 6.022141e23 sigch4 = 1.85e-17 ! aborption cross section of ch4 in cm-2 mol-1 ! Initial Solar flux Lyman alpha at Earth (1AU) flym_sol_earth=pfluxuv*1.e15 ! ph.m-2.s-1 -> about 5e+11 ph.cm-2.s-1 ! Initial Solar flux Lyman alpha on Pluto flym_sol_pluto=flym_sol_earth/(pdist_sol*pdist_sol) ! ph.m-2.s-1 ! option decrease by 12% the initial Solar Lyman alpha flux to account for Interplanetary H absorption: ! Cf Fig 3 Gladstone et al 2014 : Lyalpha at Pluto flym_sol_pluto=flym_sol_pluto*0.878 ! Time constant of conversion in aerosol [second] ! To be explored: 1.E5 - 1.E9 tcon= tcon_ch4 ! Parameter of conversion precurseur to haze kch4=k_ch4 ncratio=ncratio_ch4 ! boost for haze considering nitrogen contribution ! ratio n/c : =0.25 if there is 1N for 3C IF (firstcall) then write(*,*) 'hazecloud: haze parameters:' write(*,*) 'tcon, kch4, ncratio = ' , tcon, kch4, ncratio firstcall=.false. ENDIF ! note: mu0 = cos(solar zenith angle) ! max(mu0) at lat = declin !----------------------------------------------------------------------- ! Total IPM Lyman alpha flux at top of atmosphere ! Top of atm. Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 ! fit Fig. 4 in Randall Gladstone et al. (2015) ! Integration over semi sphere of Gladstone data. ! 145 R averaged over the sky ! 72.5 R in average over the visible hemisphere IF (ngrid.eq.1.and..not.global1d) THEN mu_sol(1) = mu0(1) mu_ipm(1) = max(mu_sol(1), 0.5) flym_ipm(1)=mu0(1)*72.5e10 ELSE IF (ngrid.eq.1.and.global1d) THEN mu_sol(1) = mu0(1) ! Full visible disk mu_ipm(1) = mu0(1)/2. ! Full sphere (day+night) flym_ipm(1)=mu_ipm(1)*145.e10 ! Proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) ? ! In theory, IPM flux remains relatively constant further than 15 AU flym_ipm(1)=flym_ipm(1)*pfluxuv/5. ELSE IF (.not.diurnal) THEN mu_sol(:) = max(mu0(:),0.) mu_ipm(:) = max(mu_sol(:), 0.5) flym_ipm(:)= mu_sol(:)*72.5e10 ELSE ! case with full fit to Gladstone et al. results ! 1) get longitude/latitude (in radian) of anti-subsolar point (max de mu0 - 180) longit=-(zday-int(zday))*2.*pi IF (longit.LE.-pi) THEN longit=longit+2.*pi ENDIF latit=-declin ! rad ! 2) Define input parameter for the fit valmin=48.74e10 ! minimum value of ipm flux in ph/m2/s valmax=115.014e10 ! maximal value of ipm flux in ph/m2/s valmin_dl=74.5e10 ! daylight minimum value puis=3.5 ! 3) Loop for each location and fit DO ig=1,ngrid ! calculation of cosinus of incident angle for IPM flux mu_sol(ig) = max(mu0(ig),0.) mu_ipm(ig) = max(mu_sol(ig), 0.5) IF (mu_sol(ig).LE.0.) THEN ! Nighttime ! Distance to subsolar point dist=acos(sin(latitude(ig))*sin(latit)+cos(latitude(ig))* & cos(latit)*cos(longit-longitude(ig)))*180/pi IF (dist.gt.90) dist=90 flym_ipm(ig)=(valmin_dl-valmin)/(90.**puis)*(dist)**puis & +valmin ELSE ! Daytime flym_ipm(ig)= mu_sol(ig)*(valmax-valmin_dl)+valmin_dl ENDIF ! Proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) ? ! In theory, IPM flux remains relatively constant further than 15 AU ! flym_ipm(ig)=flym_ipm(ig)*pfluxuv/5. ENDDO ! nlayer ENDIF ! ngrid !----------------------------------------------------------------------- ! Methane profile ! If fixed profile of CH4 gas if (haze_ch4proffix) then if (triton) then file_path=trim(datadir)//'/gas_prop/vmr_ch4_triton.txt' else file_path=trim(datadir)//'/gas_prop/vmr_ch4.txt' endif open(115,file=file_path,form='formatted') do ifine=1,Nfine read(115,*) levdat(ifine), vmrdat(ifine) enddo close(115) endif ! Initialization: vmrch4(:,:)=0. ! Diagnostics (temporary) fluxlym_sol_tot(:)=flym_sol_pluto*mu_sol(:) fluxlym_ipm_tot(:)=flym_ipm(:) call write_output("fluxlym_sol","solar lyman alpha flux","",fluxlym_sol_tot) call write_output("fluxlym_ipm","IPM lyman alpha flux","",fluxlym_ipm_tot) !----------------------------------------------------------------------- ! Main Loop DO ig=1,ngrid !! Initializations zq_ch4(ig,:)=0. zq_prec(ig,:)=0. zdqconv_prec(ig,:)=0. zdqhaze(ig,:,:)=0. zdqphot_prec(ig,:)=0. zdqphot_ch4(ig,:)=0. tauch4(:)=0. gradflux(:)=0. fluxlym_sol(1:nlayer)=0. fluxlym_sol(nlayer+1)=flym_sol_pluto*mu_sol(ig) ! fluxlym_ipm(nlayer+1)= flym_ipm(ig) !! Interpolate CH4 MMR on the model vertical grid if (haze_ch4proffix) then CALL interp_line(levdat,vmrdat,Nfine, & zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer) endif DO l=nlayer,1,-1 !! Actualisation tracer ch4 and prec_haze if (haze_ch4proffix) then zq_ch4(ig,l)=vmrch4(ig,l)*mmol(igcm_ch4_gas) / & (100.*mmol(igcm_n2)) else zq_ch4(ig,l)=pq(ig,l,igcm_ch4_gas) + & pdq(ig,l,igcm_ch4_gas)*ptimestep endif zq_prec(ig,l)=pq(ig,l,igcm_prec_haze)+ & pdq(ig,l,igcm_prec_haze)*ptimestep !! Sanity check if (zq_ch4(ig,l).lt.0.) then zq_ch4(ig,l)=0. endif if (zq_prec(ig,l).lt.0.) then zq_prec(ig,l)=0. endif !! Calculation of CH4 optical depth in Lyman alpha for each layer l tauch4(l)=sigch4*1.e-4*avogadro* & (zq_ch4(ig,l)/(mmol(igcm_ch4_gas)*1.e-3))* & (pplev(ig,l)-pplev(ig,l+1))/g !! Calculation of Flux in each layer l if (mu_sol(ig).gt.1.e-6) then fluxlym_sol(l)=fluxlym_sol(l+1)*exp(-tauch4(l)/mu_sol(ig)) endif fluxlym_ipm(l)=fluxlym_ipm(l+1)*exp(-tauch4(l)/mu_ipm(ig)) gradflux(l)=fluxlym_sol(l+1)-fluxlym_sol(l) & + fluxlym_ipm(l+1)-fluxlym_ipm(l) !! tendancy on ch4 !! considering 1 photon destroys 1 ch4 by photolysis zdqphot_ch4(ig,l)=-mmol(igcm_ch4_gas)*1.e-3/avogadro* & gradflux(l)*g/(pplev(ig,l)-pplev(ig,l+1)) !! tendency of precursors created by photolysis of ch4 zdqphot_prec(ig,l)=-zdqphot_ch4(ig,l) !! update of precursors mass: zq_prec zq_prec(ig,l)=zq_prec(ig,l)+ & zdqphot_prec(ig,l)*ptimestep !! Conversion of prec_haze into haze !! controlled by the time constant tcon !! New tendancy for prec_haze zdqconv_prec(ig,l) = -zq_prec(ig,l)* & (1-exp(-ptimestep/tcon))/ptimestep ENDDO ! nlayer !! Final tendancy for prec_haze, haze and CH4 DO iq=1,nq tname=noms(iq) if (tname(1:4).eq."haze") then zdqhaze(ig,:,iq) = -zdqconv_prec(ig,:)* & kch4*(1.+ncratio*14./12.) else if (noms(iq).eq."prec_haze") then zdqhaze(ig,:,igcm_prec_haze)= zdqphot_prec(ig,:) + & zdqconv_prec(ig,:) else if (noms(iq).eq."ch4_gas".and.(.not.hazeconservch4).and.(.not.haze_ch4proffix)) then zdqhaze(ig,:,igcm_ch4_gas)= zdqphot_ch4(ig,:) endif ENDDO ENDDO ! ngrid end subroutine hazecloud