source: trunk/LMDZ.PLUTO.old/libf/phypluto/hazecloud.F90 @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

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