| 1 | subroutine hazecloud(ngrid,nlayer,nq,ptimestep, & |
|---|
| 2 | pplay,pplev,pq,pdq,pdist_sol,mu0,pfluxuv,zdqhaze, & |
|---|
| 3 | zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin) |
|---|
| 4 | |
|---|
| 5 | ! zdqphot_ch4,zdqconv_prec,declin,zdqhaze_col) |
|---|
| 6 | |
|---|
| 7 | use radinc_h, only : naerkind |
|---|
| 8 | use comgeomfi_h |
|---|
| 9 | use comcstfi_mod, only: pi, g |
|---|
| 10 | use tracer_h, only: igcm_haze, igcm_ch4_gas, igcm_prec_haze, noms, mmol |
|---|
| 11 | use geometry_mod, only: longitude, latitude ! in radians |
|---|
| 12 | use callkeys_mod, only: hazeconservch4 |
|---|
| 13 | |
|---|
| 14 | implicit none |
|---|
| 15 | |
|---|
| 16 | !================================================================== |
|---|
| 17 | ! Purpose |
|---|
| 18 | ! ------- |
|---|
| 19 | ! Production of haze in the atmosphere |
|---|
| 20 | ! |
|---|
| 21 | ! Inputs |
|---|
| 22 | ! ------ |
|---|
| 23 | ! ngrid Number of vertical columns |
|---|
| 24 | ! nlayer Number of layers |
|---|
| 25 | ! pplay(ngrid,nlayer) Pressure layers |
|---|
| 26 | ! pplev(ngrid,nlayer+1) Pressure levels |
|---|
| 27 | ! |
|---|
| 28 | ! Outputs |
|---|
| 29 | ! ------- |
|---|
| 30 | ! |
|---|
| 31 | ! Both |
|---|
| 32 | ! ---- |
|---|
| 33 | ! |
|---|
| 34 | ! Authors |
|---|
| 35 | ! ------- |
|---|
| 36 | ! Tanguy Bertrand and Francois Forget (2014) |
|---|
| 37 | ! |
|---|
| 38 | !================================================================== |
|---|
| 39 | |
|---|
| 40 | |
|---|
| 41 | !----------------------------------------------------------------------- |
|---|
| 42 | ! Arguments |
|---|
| 43 | |
|---|
| 44 | INTEGER ngrid, nlayer, nq |
|---|
| 45 | ! REAL lati(ngrid) |
|---|
| 46 | ! REAL long(ngrid) |
|---|
| 47 | ! REAL declin |
|---|
| 48 | LOGICAL firstcall |
|---|
| 49 | SAVE firstcall |
|---|
| 50 | DATA firstcall/.true./ |
|---|
| 51 | REAL ptimestep |
|---|
| 52 | REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1) |
|---|
| 53 | REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) |
|---|
| 54 | REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) |
|---|
| 55 | ! REAL,INTENT(IN) :: mmol(nq) |
|---|
| 56 | REAL,INTENT(IN) :: pdist_sol ! distance SUN-pluto in AU |
|---|
| 57 | REAL,INTENT(IN) :: pfluxuv ! Lyman alpha flux at specific Ls (ph/cm/s) |
|---|
| 58 | REAL,INTENT(IN) :: mu0(ngrid) ! cosinus of solar incident flux |
|---|
| 59 | REAL,INTENT(IN) :: declin ! distance SUN-pluto in AU |
|---|
| 60 | REAL,INTENT(OUT) :: zdqhaze(ngrid,nlayer,nq) ! Final tendancy |
|---|
| 61 | ! REAL,INTENT(OUT) :: zdqhaze_col(ngrid) ! Final tendancy haze prod |
|---|
| 62 | !----------------------------------------------------------------------- |
|---|
| 63 | ! Local variables |
|---|
| 64 | |
|---|
| 65 | INTEGER l,ig,iq |
|---|
| 66 | REAL zq_ch4(ngrid,nlayer) |
|---|
| 67 | REAL zq_prec(ngrid,nlayer) |
|---|
| 68 | REAL tauch4(nlayer) |
|---|
| 69 | REAL sigch4 ! aborption cross section ch4 cm-2 mol-1 |
|---|
| 70 | REAL flym_sol_earth ! initial flux Earth ph.m-2.s-1 |
|---|
| 71 | REAL flym_sol_pluto ! initial Lyman alpha SOLAR flux Pluto ph.m-2.s-1 |
|---|
| 72 | REAL flym_ipm(ngrid) ! top of atm. Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 |
|---|
| 73 | REAL fluxlym_sol(nlayer+1) ! Local Solar flux Lyman alpha ph.m-2.s-1 |
|---|
| 74 | REAL fluxlym_ipm(nlayer+1) ! Local IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 |
|---|
| 75 | REAL mu_ipm(ngrid) ! local Mean incident flux for IPM Lyman alpha photons |
|---|
| 76 | REAL mu_sol(ngrid) ! local Mean incident flux for solar lyman alpha photons |
|---|
| 77 | |
|---|
| 78 | REAL gradflux(nlayer) ! gradient flux Lyman alpha ph.m-2.s-1 |
|---|
| 79 | REAL zdqphot_ch4(ngrid,nlayer) ! tendancy due to photolysis of ch4 |
|---|
| 80 | REAL zdqphot_prec(ngrid,nlayer) ! tendancy due to photolysis of ch4 |
|---|
| 81 | REAL zdqconv_prec(ngrid,nlayer) ! tendancy due to conversion of |
|---|
| 82 | ! prec_haze into haze |
|---|
| 83 | REAL kch4 ! fraction of Carbon from ch4 directly dissociated |
|---|
| 84 | ! by prec_haze |
|---|
| 85 | CHARACTER(len=10) :: tname |
|---|
| 86 | REAL ncratio ! ration Nitrogen/Carbon observed in tholins |
|---|
| 87 | REAL tcon ! Time constant: conversion in aerosol |
|---|
| 88 | REAL avogadro ! avogadro constant |
|---|
| 89 | |
|---|
| 90 | REAL longit |
|---|
| 91 | REAL latit |
|---|
| 92 | REAL valmin |
|---|
| 93 | REAL valmax |
|---|
| 94 | REAL valmin_dl |
|---|
| 95 | REAL puis |
|---|
| 96 | REAL dist |
|---|
| 97 | REAL long2 |
|---|
| 98 | !----------------------------------------------------------------------- |
|---|
| 99 | |
|---|
| 100 | !---------------- INPUT ------------------------------------------------ |
|---|
| 101 | |
|---|
| 102 | avogadro = 6.022141e23 |
|---|
| 103 | sigch4 = 1.85e-17 ! aborption cross section ch4 cm-2 mol-1 |
|---|
| 104 | !! Initial Solar flux Lyman alpha on Earth (1AU) |
|---|
| 105 | flym_sol_earth=pfluxuv*1.e15 ! ph.m-2.s-1 -> 5e+11 ph.cm-2.s-1 |
|---|
| 106 | !! Initial Solar flux Lyman alpha on Pluto |
|---|
| 107 | flym_sol_pluto=flym_sol_earth/(pdist_sol*pdist_sol) ! ph.m-2.s-1 |
|---|
| 108 | ! option decrease by 12% the initial IPM flux to account for Interplanetary H absorption: |
|---|
| 109 | ! Fig 3 Gladstone et al 2014 : Lyalpha at Pluto |
|---|
| 110 | flym_sol_pluto=flym_sol_pluto*0.878 |
|---|
| 111 | |
|---|
| 112 | !!!! Top of atm. Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1 |
|---|
| 113 | ! fit Fig. 4 in Randall Gladstone et al. (2015) |
|---|
| 114 | ! -- New version : Integration over semi sphere of Gladstone data. |
|---|
| 115 | ! Fit of results |
|---|
| 116 | |
|---|
| 117 | IF (ngrid.eq.1) THEN |
|---|
| 118 | flym_ipm(1)=75.e10 |
|---|
| 119 | mu_ipm(1) = 0.5 !max(mu0(1), 0.5) |
|---|
| 120 | mu_sol(1)=0.25 |
|---|
| 121 | ELSE |
|---|
| 122 | |
|---|
| 123 | ! 1) get longitude/latitude (rad) of anti-subsolar point (max de mu0 - 180) |
|---|
| 124 | longit=longitude((MAXLOC(mu0,DIM=1,MASK=mu0.GT.0.9))) ! rad |
|---|
| 125 | IF (longit.GE.0) THEN |
|---|
| 126 | longit=longit-pi |
|---|
| 127 | ELSE |
|---|
| 128 | longit=longit+pi |
|---|
| 129 | ENDIF |
|---|
| 130 | latit=-declin ! rad |
|---|
| 131 | |
|---|
| 132 | ! 2) Define input parameter for the fit |
|---|
| 133 | valmin=48.74e10 ! minimum value of ipm flux in ph/m2/s |
|---|
| 134 | valmax=115.014e10 ! maximal value of ipm flux in ph/m2/s |
|---|
| 135 | valmin_dl=74.5e10 ! daylight minimum value |
|---|
| 136 | puis=3.5 |
|---|
| 137 | |
|---|
| 138 | ! 3) Loop for each location and fit |
|---|
| 139 | DO ig=1,ngrid |
|---|
| 140 | ! calculation of cosinus of incident angle for IPM flux |
|---|
| 141 | mu_sol(ig) = mu0(ig) |
|---|
| 142 | mu_ipm(ig) = max(mu_sol(ig), 0.5) |
|---|
| 143 | IF (mu0(ig).LT.1.e-4) THEN |
|---|
| 144 | dist=acos(sin(latitude(ig))*sin(latit)+cos(latitude(ig))* & |
|---|
| 145 | cos(latit)*cos(longit-longitude(ig)))*180/pi |
|---|
| 146 | IF (dist.gt.90) dist=90 |
|---|
| 147 | flym_ipm(ig)=(valmin_dl-valmin)/(90.**puis)*(dist)**puis & |
|---|
| 148 | +valmin |
|---|
| 149 | ELSE |
|---|
| 150 | flym_ipm(ig)= mu0(ig)*(valmax-valmin_dl)+valmin_dl |
|---|
| 151 | ENDIF |
|---|
| 152 | ! proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) |
|---|
| 153 | flym_ipm(ig)=flym_ipm(ig)*pfluxuv/5. |
|---|
| 154 | ENDDO ! nlayer |
|---|
| 155 | |
|---|
| 156 | |
|---|
| 157 | ENDIF ! ngrid |
|---|
| 158 | !--- |
|---|
| 159 | |
|---|
| 160 | !! Time constant of conversion in aerosol to be explore |
|---|
| 161 | tcon= 1.e7 ! 153 ! 1.E7 !second |
|---|
| 162 | ! tcon= 1. ! 153 ! 1.E7 !second |
|---|
| 163 | !! Parameter of conversion precurseur to haze |
|---|
| 164 | kch4=1. |
|---|
| 165 | ncratio=0.5 ! boost for haze considering nitrogen contribution |
|---|
| 166 | ! ratio n/c : =0.25 if there is 1N for 3C |
|---|
| 167 | IF (firstcall) then |
|---|
| 168 | write(*,*) 'hazecloud: haze parameters:' |
|---|
| 169 | write(*,*) 'tcon, kch4, ncratio = ' , tcon, kch4, ncratio |
|---|
| 170 | firstcall=.false. |
|---|
| 171 | ENDIF |
|---|
| 172 | ! note: mu0 = cos(solar zenith angle) |
|---|
| 173 | ! max(mu0) = declin |
|---|
| 174 | |
|---|
| 175 | !!---------------------------------------------------------- |
|---|
| 176 | !!---------------------------------------------------------- |
|---|
| 177 | |
|---|
| 178 | DO ig=1,ngrid |
|---|
| 179 | |
|---|
| 180 | zq_ch4(ig,:)=0. |
|---|
| 181 | zq_prec(ig,:)=0. |
|---|
| 182 | zdqconv_prec(ig,:)=0. |
|---|
| 183 | zdqhaze(ig,:,:)=0. |
|---|
| 184 | zdqphot_prec(ig,:)=0. |
|---|
| 185 | zdqphot_ch4(ig,:)=0. |
|---|
| 186 | tauch4(:)=0. |
|---|
| 187 | gradflux(:)=0. |
|---|
| 188 | fluxlym_sol(1:nlayer)=0. |
|---|
| 189 | fluxlym_sol(nlayer+1)=flym_sol_pluto*mu_sol(ig) ! |
|---|
| 190 | fluxlym_ipm(nlayer+1)= flym_ipm(ig) |
|---|
| 191 | |
|---|
| 192 | DO l=nlayer,1,-1 |
|---|
| 193 | !! Actualisation tracer ch4 and prec_haze |
|---|
| 194 | |
|---|
| 195 | !IF (ngrid.eq.1) THEN |
|---|
| 196 | !! option zq_ch4 = cte |
|---|
| 197 | ! zq_ch4(ig,l)=0.01*0.5*16./28. ! Temporaire |
|---|
| 198 | !ELSE |
|---|
| 199 | zq_ch4(ig,l)=pq(ig,l,igcm_ch4_gas)+ & |
|---|
| 200 | pdq(ig,l,igcm_ch4_gas)*ptimestep |
|---|
| 201 | !ENDIF |
|---|
| 202 | if (zq_ch4(ig,l).lt.0.) then |
|---|
| 203 | zq_ch4(ig,l)=0. |
|---|
| 204 | endif |
|---|
| 205 | |
|---|
| 206 | !! option zq_ch4 = cte |
|---|
| 207 | ! zq_ch4(ig,l)=0.1*16./28. ! Temporaire |
|---|
| 208 | |
|---|
| 209 | zq_prec(ig,l)=pq(ig,l,igcm_prec_haze)+ & |
|---|
| 210 | pdq(ig,l,igcm_prec_haze)*ptimestep |
|---|
| 211 | |
|---|
| 212 | !! Calculation optical depth ch4 in Lyman alpha for each layer l |
|---|
| 213 | tauch4(l)=sigch4*1.e-4*avogadro* & |
|---|
| 214 | (zq_ch4(ig,l)/(mmol(igcm_ch4_gas)*1.e-3))* & |
|---|
| 215 | (pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 216 | !! Calculation of Flux in each layer l |
|---|
| 217 | if (mu_sol(ig).gt.1.e-6) then |
|---|
| 218 | fluxlym_sol(l)=fluxlym_sol(l+1)*exp(-tauch4(l)/mu_sol(ig)) |
|---|
| 219 | endif |
|---|
| 220 | |
|---|
| 221 | fluxlym_ipm(l)=fluxlym_ipm(l+1)*exp(-tauch4(l)/mu_ipm(ig)) |
|---|
| 222 | |
|---|
| 223 | gradflux(l)=fluxlym_sol(l+1)-fluxlym_sol(l) & |
|---|
| 224 | + fluxlym_ipm(l+1)-fluxlym_ipm(l) |
|---|
| 225 | !! tendancy on ch4 |
|---|
| 226 | !! considering 1 photon destroys 1 ch4 by photolysis |
|---|
| 227 | zdqphot_ch4(ig,l)=-mmol(igcm_ch4_gas)*1.e-3/avogadro* & |
|---|
| 228 | gradflux(l)*g/(pplev(ig,l)-pplev(ig,l+1)) |
|---|
| 229 | |
|---|
| 230 | !! tendency of prec created by photolysis of ch4 |
|---|
| 231 | zdqphot_prec(ig,l)=-zdqphot_ch4(ig,l) |
|---|
| 232 | |
|---|
| 233 | !! update precurseur zq_prec |
|---|
| 234 | zq_prec(ig,l)=zq_prec(ig,l)+ & |
|---|
| 235 | zdqphot_prec(ig,l)*ptimestep |
|---|
| 236 | |
|---|
| 237 | !! Conversion of prec_haze into haze |
|---|
| 238 | !! controlled by the time constant tcon |
|---|
| 239 | !! New tendancy for prec_haze |
|---|
| 240 | zdqconv_prec(ig,l) = -zq_prec(ig,l)* & |
|---|
| 241 | (1-exp(-ptimestep/tcon))/ptimestep |
|---|
| 242 | |
|---|
| 243 | ENDDO ! nlayer |
|---|
| 244 | |
|---|
| 245 | !! Final tendancy for prec_haze and haze |
|---|
| 246 | |
|---|
| 247 | DO iq=1,nq |
|---|
| 248 | tname=noms(iq) |
|---|
| 249 | !print*, 'TB17 tname=',tname,tname(1:4) |
|---|
| 250 | if (tname(1:4).eq."haze") then |
|---|
| 251 | zdqhaze(ig,:,iq) = -zdqconv_prec(ig,:)* & |
|---|
| 252 | kch4*(1.+ncratio*14./12.) |
|---|
| 253 | else if (noms(iq).eq."prec_haze") then |
|---|
| 254 | zdqhaze(ig,:,igcm_prec_haze)= zdqphot_prec(ig,:) + & |
|---|
| 255 | zdqconv_prec(ig,:) |
|---|
| 256 | |
|---|
| 257 | else if (noms(iq).eq."ch4_gas".and.(.not.hazeconservch4)) then |
|---|
| 258 | zdqhaze(ig,:,igcm_ch4_gas)= zdqphot_ch4(ig,:) |
|---|
| 259 | endif |
|---|
| 260 | ENDDO |
|---|
| 261 | |
|---|
| 262 | ENDDO ! ngrid |
|---|
| 263 | |
|---|
| 264 | |
|---|
| 265 | !! tendency kg/m2/s for haze column: |
|---|
| 266 | ! zdqhaze_col(:)=0. |
|---|
| 267 | ! DO ig=1,ngrid |
|---|
| 268 | ! DO l=1,nlayer |
|---|
| 269 | ! zdqhaze_col(ig)=zdqhaze_col(ig)+zdqhaze(ig,l,igcm_haze)* & |
|---|
| 270 | ! (pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 271 | ! ENDDO |
|---|
| 272 | ! ENDDO |
|---|
| 273 | |
|---|
| 274 | end subroutine hazecloud |
|---|
| 275 | |
|---|
