| 1 | subroutine prodhaze(ngrid,nlayer,nq,ptimestep, & |
|---|
| 2 | pplev,pq,pdq,pdist_sol,mu0,declin,zdqfasthaze,zdqsfasthaze,gradflux, & |
|---|
| 3 | flym_bot,flym_sol_bot,flym_ipm_bot,flym_sol,flym_ipm) |
|---|
| 4 | |
|---|
| 5 | use comgeomfi_h |
|---|
| 6 | use comcstfi_mod, only: pi, g |
|---|
| 7 | use tracer_h, only: igcm_ch4_gas, igcm_prec_haze, noms, mmol |
|---|
| 8 | use geometry_mod, only: longitude, latitude ! in radians |
|---|
| 9 | use callkeys_mod, only: hazeconservch4, diurnal |
|---|
| 10 | implicit none |
|---|
| 11 | |
|---|
| 12 | !================================================================== |
|---|
| 13 | ! Purpose |
|---|
| 14 | ! ------- |
|---|
| 15 | ! Fast production of haze in the atmosphere, direct accumulation to surface |
|---|
| 16 | ! |
|---|
| 17 | ! Authors |
|---|
| 18 | ! ------- |
|---|
| 19 | ! Tanguy Bertrand and Francois Forget (2014) |
|---|
| 20 | ! |
|---|
| 21 | !================================================================== |
|---|
| 22 | |
|---|
| 23 | !----------------------------------------------------------------------- |
|---|
| 24 | ! Arguments |
|---|
| 25 | |
|---|
| 26 | INTEGER ngrid, nlayer,nq |
|---|
| 27 | LOGICAL firstcall |
|---|
| 28 | SAVE firstcall |
|---|
| 29 | DATA firstcall/.true./ |
|---|
| 30 | REAL ptimestep |
|---|
| 31 | REAL pplev(ngrid,nlayer+1) |
|---|
| 32 | REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) |
|---|
| 33 | REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) |
|---|
| 34 | REAL,INTENT(IN) :: pdist_sol ! distance SUN-pluto in AU |
|---|
| 35 | REAL,INTENT(IN) :: mu0(ngrid) ! cosinus of solar incident flux |
|---|
| 36 | REAL,INTENT(IN) :: declin ! declination |
|---|
| 37 | REAL,INTENT(OUT) :: zdqfasthaze(ngrid,nq) ! Final tendancy ch4 |
|---|
| 38 | REAL,INTENT(OUT) :: zdqsfasthaze(ngrid) ! Final tendancy haze surf |
|---|
| 39 | REAL, INTENT(OUT) :: gradflux(ngrid) ! gradient flux Lyman alpha ph.m-2.s-1 |
|---|
| 40 | REAL, INTENT(OUT) :: flym_bot(ngrid) ! flux Lyman alpha bottom ph.m-2.s-1 |
|---|
| 41 | REAL, INTENT(OUT) :: flym_sol_bot(ngrid) ! Solar flux Lyman alpha ph.m-2.s-1 reaching the surface |
|---|
| 42 | REAL, INTENT(OUT) :: flym_ipm_bot(ngrid) ! IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 reaching the surface |
|---|
| 43 | REAL, INTENT(OUT) :: flym_sol(ngrid) ! Incident Solar flux Lyman alpha ph.m-2.s-1 |
|---|
| 44 | REAL, INTENT(OUT) :: flym_ipm(ngrid) ! Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 |
|---|
| 45 | !----------------------------------------------------------------------- |
|---|
| 46 | ! Local variables |
|---|
| 47 | |
|---|
| 48 | INTEGER l,ig,iq |
|---|
| 49 | REAL zq_ch4(ngrid) |
|---|
| 50 | REAL zq_prec(ngrid) |
|---|
| 51 | REAL tauch4 |
|---|
| 52 | REAL sigch4 ! aborption cross section ch4 cm-2 mol-1 |
|---|
| 53 | REAL flym_sol_earth ! initial flux Earth ph.m-2.s-1 |
|---|
| 54 | REAL flym_sol_pluto ! initial Lyman alpha SOLAR flux Pluto ph.m-2.s-1 |
|---|
| 55 | REAL mu_ipm(ngrid) ! local Mean incident flux for IPM Lyman alpha photons |
|---|
| 56 | |
|---|
| 57 | REAL kch4 ! fraction of Carbon from ch4 directly dissociated by prec_haze |
|---|
| 58 | REAL ncratio ! ration Nitrogen/Carbon observed in tholins |
|---|
| 59 | REAL avogadro ! avogadro constant |
|---|
| 60 | |
|---|
| 61 | CHARACTER(len=10) :: tname |
|---|
| 62 | REAL pagg,pform |
|---|
| 63 | REAL tconv |
|---|
| 64 | REAL longit |
|---|
| 65 | REAL latit |
|---|
| 66 | REAL valmin |
|---|
| 67 | REAL valmax |
|---|
| 68 | REAL valmin_dl |
|---|
| 69 | REAL puis |
|---|
| 70 | REAL dist |
|---|
| 71 | REAL zdqprec(ngrid) |
|---|
| 72 | REAL zdqch4(ngrid) |
|---|
| 73 | REAL zdqconv_prec(ngrid) |
|---|
| 74 | !----------------------------------------------------------------------- |
|---|
| 75 | |
|---|
| 76 | !******** INPUT |
|---|
| 77 | avogadro = 6.022141e23 |
|---|
| 78 | sigch4 = 1.85e-17 ! aborption cross section ch4 cm-2 mol-1 |
|---|
| 79 | !! Initial Solar flux Lyman alpha on Earth (1AU) |
|---|
| 80 | flym_sol_earth=5.e15 ! ph.m-2.s-1 --> 5e+11 ph.cm-2.s-1 |
|---|
| 81 | !! Parameter of conversion precurseur to haze : here accelaeration/deceleration of haze producition |
|---|
| 82 | kch4=1. |
|---|
| 83 | ncratio=0.5 ! boost for haze considering nitrogen contribution |
|---|
| 84 | ! ratio n/c : =0.25 if there is 1N for 3C |
|---|
| 85 | pagg=0.05 ! Minimum pressure required for aggregation |
|---|
| 86 | pform=0.006 ! Minimum pressure required for haze formation |
|---|
| 87 | tconv=1.e7 ! (s) time needed to convert gas precursors into aerosol |
|---|
| 88 | !******** Incident Solar flux |
|---|
| 89 | !! Initial Solar flux Lyman alpha on Pluto |
|---|
| 90 | flym_sol_pluto=0.25*flym_sol_earth/(pdist_sol*pdist_sol) ! ph.m-2.s-1 |
|---|
| 91 | ! Decreasing by 12% the initial IPM flux to account for Interplanetary H absorption: |
|---|
| 92 | ! Fig 3 Gladstone et al 2014 : Lyalpha at Pluto |
|---|
| 93 | flym_sol_pluto=flym_sol_pluto*(1.-0.12*pdist_sol/32.91) |
|---|
| 94 | |
|---|
| 95 | !******** Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 |
|---|
| 96 | ! fit Fig. 4 in Randall Gladstone et al. (2015) |
|---|
| 97 | ! -- New version : Integration over semi sphere of Gladstone data. |
|---|
| 98 | ! Fit of results |
|---|
| 99 | |
|---|
| 100 | if (diurnal) then |
|---|
| 101 | ! 1) get longitude/latitude (rad) of anti-subsolar point (max de mu0 - 180) |
|---|
| 102 | longit=longitude((MAXLOC(mu0,DIM=1,MASK=mu0.GT.0.9))) ! rad |
|---|
| 103 | IF (longit.GE.0) THEN |
|---|
| 104 | longit=longit-pi |
|---|
| 105 | ELSE |
|---|
| 106 | longit=longit+pi |
|---|
| 107 | ENDIF |
|---|
| 108 | latit=-declin ! rad |
|---|
| 109 | |
|---|
| 110 | ! 2) Define input parameter for the fit |
|---|
| 111 | valmin=48.74e10 ! minimum value of ipm flux in ph/m2/s |
|---|
| 112 | valmax=115.014e10 ! maximal value of ipm flux in ph/m2/s |
|---|
| 113 | valmin_dl=74.5e10 ! daylight minimum value |
|---|
| 114 | puis=3.5 |
|---|
| 115 | |
|---|
| 116 | ! 3) Loop for each location and fit |
|---|
| 117 | DO ig=1,ngrid |
|---|
| 118 | !!! Solar flux |
|---|
| 119 | flym_sol(ig)=flym_sol_pluto*mu0(ig) |
|---|
| 120 | |
|---|
| 121 | !!! IPM flux |
|---|
| 122 | ! calculation of cosinus of incident angle for IPM flux |
|---|
| 123 | mu_ipm(ig) = max(mu0(ig), 0.5) |
|---|
| 124 | IF (mu0(ig).LT.1.e-4) THEN |
|---|
| 125 | dist=acos(sin(latitude(ig))*sin(latit)+cos(latitude(ig))* & |
|---|
| 126 | cos(latit)*cos(longit-longitude(ig)))*180/pi |
|---|
| 127 | IF (dist.gt.90) dist=90 |
|---|
| 128 | flym_ipm(ig)=(valmin_dl-valmin)/(90.**puis)*(dist)**puis & |
|---|
| 129 | +valmin |
|---|
| 130 | ELSE |
|---|
| 131 | flym_ipm(ig)= mu0(ig)*(valmax-valmin_dl)+valmin_dl |
|---|
| 132 | |
|---|
| 133 | !! IPM flux must be weighted by solar flux |
|---|
| 134 | flym_ipm(ig)=flym_ipm(ig)*flym_sol_pluto/1.0156e12 ! Value in 2015 for flym_sol_pluto |
|---|
| 135 | ENDIF |
|---|
| 136 | ENDDO ! ngrid |
|---|
| 137 | |
|---|
| 138 | ELSE ! diurnal=false |
|---|
| 139 | |
|---|
| 140 | DO ig=1,ngrid |
|---|
| 141 | flym_sol(ig)=flym_sol_pluto*mu0(ig) |
|---|
| 142 | flym_ipm(ig)=flym_sol(ig)*1.15 |
|---|
| 143 | mu_ipm(ig) = max(mu0(ig), 0.5) |
|---|
| 144 | ENDDO |
|---|
| 145 | |
|---|
| 146 | ENDIF ! diurnal |
|---|
| 147 | |
|---|
| 148 | IF (firstcall) then |
|---|
| 149 | write(*,*) 'Fast haze parameters:' |
|---|
| 150 | write(*,*) 'kch4, ncratio = ',kch4,ncratio |
|---|
| 151 | firstcall=.false. |
|---|
| 152 | ENDIF |
|---|
| 153 | |
|---|
| 154 | ! note: mu0 = cos(solar zenith angle) |
|---|
| 155 | ! max(mu0) = declin |
|---|
| 156 | |
|---|
| 157 | !************************** |
|---|
| 158 | !******** LOOP ************ |
|---|
| 159 | !************************** |
|---|
| 160 | zdqfasthaze(:,:)=0. |
|---|
| 161 | zdqsfasthaze(:)=0. |
|---|
| 162 | zdqch4(:)=0. |
|---|
| 163 | zdqprec(:)=0. |
|---|
| 164 | DO ig=1,ngrid |
|---|
| 165 | |
|---|
| 166 | !! Update of tracer ch4 and prec_haze |
|---|
| 167 | zq_ch4(ig)=pq(ig,1,igcm_ch4_gas)+pdq(ig,1,igcm_ch4_gas)*ptimestep |
|---|
| 168 | zq_prec(ig)=pq(ig,1,igcm_prec_haze)+pdq(ig,1,igcm_prec_haze)*ptimestep |
|---|
| 169 | !! Security as we loose CH4 molecules |
|---|
| 170 | if (zq_ch4(ig).lt.0.) then |
|---|
| 171 | zq_ch4(ig)=0. |
|---|
| 172 | endif |
|---|
| 173 | |
|---|
| 174 | !! Calculation of CH4 optical depth in Lyman alpha for whole atm |
|---|
| 175 | tauch4=sigch4*1.e-4*avogadro*(zq_ch4(ig)/(mmol(igcm_ch4_gas)*1.e-3))*(pplev(ig,1))/g |
|---|
| 176 | |
|---|
| 177 | !! Calculation of Flux reaching the surface |
|---|
| 178 | if (mu0(ig).gt.1.e-6) then |
|---|
| 179 | flym_sol_bot(ig)=flym_sol(ig)*exp(-tauch4/(mu0(ig))) |
|---|
| 180 | endif |
|---|
| 181 | flym_ipm_bot(ig)=flym_ipm(ig)*exp(-tauch4/mu_ipm(ig)) |
|---|
| 182 | |
|---|
| 183 | !! Total flux reaching the surface and flux absorbed by CH4 |
|---|
| 184 | gradflux(ig)=flym_sol(ig)-flym_sol_bot(ig) + flym_ipm(ig)-flym_ipm_bot(ig) |
|---|
| 185 | flym_bot(ig)=flym_ipm_bot(ig)+flym_sol_bot(ig) |
|---|
| 186 | |
|---|
| 187 | !! Tendancy on ch4 considering 1 photon destroys 1 ch4 by photolysis |
|---|
| 188 | zdqch4(ig)=-mmol(igcm_ch4_gas)*1.e-3/avogadro*gradflux(ig)*g/(pplev(ig,1)) |
|---|
| 189 | !! Tendency of precursors created |
|---|
| 190 | zdqprec(ig)=-zdqch4(ig) |
|---|
| 191 | |
|---|
| 192 | !! update precurseur zq_prec |
|---|
| 193 | zq_prec(ig)=zq_prec(ig)+zdqprec(ig)*ptimestep |
|---|
| 194 | |
|---|
| 195 | !! If pressure > pagg, aggregation is possible |
|---|
| 196 | !! If pressure > pform, haze formation is possible |
|---|
| 197 | if (pplev(ig,1).gt.pform) then |
|---|
| 198 | ! New tendancy for prec_haze, turning into haze aerosols |
|---|
| 199 | zdqconv_prec(ig) = -zq_prec(ig)*(1.-exp(-ptimestep/tconv))/ptimestep |
|---|
| 200 | else |
|---|
| 201 | zdqconv_prec(ig) = 0. |
|---|
| 202 | endif |
|---|
| 203 | |
|---|
| 204 | !! Final tendencies |
|---|
| 205 | DO iq=1,nq |
|---|
| 206 | tname=noms(iq) |
|---|
| 207 | if (tname(1:4).eq."haze") then |
|---|
| 208 | ! Tendency haze in atmosphere kg/kg |
|---|
| 209 | zdqfasthaze(ig,iq) = -zdqconv_prec(ig)*kch4*(1.+ncratio*14./12.) |
|---|
| 210 | ! Tendancy for haze directly deposited at the surface kg/m2 |
|---|
| 211 | zdqsfasthaze(ig) = zdqfasthaze(ig,iq)*pplev(ig,1)/g |
|---|
| 212 | else if (noms(iq).eq."prec_haze") then |
|---|
| 213 | zdqfasthaze(ig,igcm_prec_haze)= zdqprec(ig) + zdqconv_prec(ig) |
|---|
| 214 | else if (noms(iq).eq."ch4_gas".and.(.not.hazeconservch4)) then |
|---|
| 215 | ! Tendancy for haze directly deposited at the surface kg/m2 |
|---|
| 216 | zdqfasthaze(ig,igcm_ch4_gas)= zdqch4(ig) |
|---|
| 217 | endif |
|---|
| 218 | ENDDO |
|---|
| 219 | |
|---|
| 220 | |
|---|
| 221 | |
|---|
| 222 | ENDDO ! ngrid |
|---|
| 223 | |
|---|
| 224 | end subroutine prodhaze |
|---|
| 225 | |
|---|
