source: trunk/mesoscale/LMDZ.MARS/libf_gcm/phymars/dustopacity.F @ 57

Last change on this file since 57 was 57, checked in by aslmd, 14 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)
File size: 9.1 KB
Line 
1      SUBROUTINE dustopacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,pq,
2     $    tauref,tau,aerosol)
3                                                   
4       IMPLICIT NONE
5c=======================================================================
6c   subject:
7c   --------
8c   Computing aerosol optical depth (dust opacity)
9c   In each layers
10c
11c   author: F.Forget
12c   ------
13c   update F. Montmessin (water ice scheme)
14c      and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
15c   
16c   input:
17c   -----
18c   ngrid             Number of gridpoint of horizontal grid
19c   nlayer            Number of layer
20c   nq                Number of tracer
21c   ls                Solar longitude (Ls) , radian
22c   pplay,pplev       pressure (Pa) in the middle and boundary of each layer
23c   pq                Dust mixing ratio (used if tracer =T and active=T).
24c
25c   output:
26c   -------
27c   tauref       Prescribed mean column optical depth at 700 Pa
28c   tau          Column total visible dust optical depth at each point
29c   aerosol      aerosol(ig,l,1) is the dust optical
30c                depth in layer l, grid point ig
31
32c
33c=======================================================================
34#include "dimensions.h"
35#include "dimphys.h"
36#include "callkeys.h"
37#include "comcstfi.h"
38#include "comgeomfi.h"
39#include "dimradmars.h"
40#include "yomaer.h"
41#include "tracer.h"
42#include "planete.h"
43
44c-----------------------------------------------------------------------
45c
46c    Declarations :
47c    --------------
48c
49c    Input/Output
50c    ------------
51      INTEGER ngrid,nlayer,nq
52
53      REAL ls,zday,expfactor   
54      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
55      REAL pq(ngrid,nlayer,nq)
56      REAL tauref(ngrid), tau(ngrid,naerkind)
57      REAL aerosol(ngrid,nlayer,naerkind)
58c
59c    Local variables :
60c    -----------------
61      INTEGER l,ig,iq
62      real topdust(ngridmx)
63      real zlsconst, zp
64      real taueq,tauS,tauN
65      real r0,reff,coefsize
66c
67c   local saved variables
68c   ---------------------
69
70      REAL topdust0(ngridmx)
71      SAVE topdust0
72
73      LOGICAL firstcall
74      DATA firstcall/.true./
75      SAVE firstcall
76
77
78c----------------------------------------------------------------------
79
80c     Initialisation
81c     --------------
82
83      IF (firstcall) THEN
84
85c        altitude of the top of the aerosol layer (km) at Ls=2.76rad:
86c        in the Viking year scenario
87         DO ig=1,ngrid
88             topdust0(ig)=60. -22.*SIN(lati(ig))**2
89         END DO
90         firstcall=.false.
91      END IF
92
93c  -------------------------------------------------------------
94c     1) Prescribed dust  (if tracer=F or active=F)
95c  -------------------------------------------------------------
96      IF ((.not.tracer) .or. (.not.active)) THEN
97
98c       Vertical column optical depth at 700.Pa
99c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
100        IF(iaervar.eq.1) THEN
101           do ig=1, ngridmx
102            tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste as read in starfi
103          end do
104        ELSE IF (iaervar.eq.2) THEN   ! << "Viking" Scenario>>
105
106          tauref(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1
107          do ig=2,ngrid
108            tauref(ig) = tauref(1)
109          end do
110
111        ELSE IF (iaervar.eq.3) THEN  ! << "MGS" scenario >>
112
113           taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls-4.363)))**14
114           tauS= 0.1 +(0.5-0.1)  *(cos(0.5*(ls-4.363)))**14
115           tauN = 0.1
116c          if (peri_day.eq.150) then
117c            tauS=0.1
118c            tauN=0.1 +(0.5-0.1)  *(cos(0.5*(ls+pi-4.363)))**14
119c            taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls+pi-4.363)))**14
120c           endif
121             
122           do ig=1,ngrid/2  ! Northern hemisphere
123             tauref(ig)= tauN +
124     &       (taueq-tauN)*0.5*(1+tanh((45-lati(ig)*180./pi)*6/60))
125           end do
126           do ig=ngrid/2+1, ngridmx  ! Southern hemisphere
127             tauref(ig)= tauS +
128     &       (taueq-tauS)*0.5*(1+tanh((45+lati(ig)*180./pi)*6/60))
129           end do
130
131        ELSE IF (iaervar.eq.4) THEN  ! << "TES scenario >>
132           call readtesassim(ngrid,nlayer,zday,pplev,tauref)
133
134        ELSE IF (iaervar.eq.5) THEN   ! << Escalier Scenario>>
135c         tauref(1) = 0.2
136c         if ((ls.ge.210.*pi/180.).and.(ls.le.330.*pi/180.))
137c    &                              tauref(1) = 2.5
138          tauref(1) = 2.5
139          if ((ls.ge.30.*pi/180.).and.(ls.le.150.*pi/180.))
140     &                              tauref(1) = .2
141
142          do ig=2,ngrid
143            tauref(ig) = tauref(1)
144          end do
145        ELSE
146          stop 'problem with iaervar in dustopacity.F'
147        ENDIF
148
149c       Altitude of the top of the dust layer
150c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
151        zlsconst=SIN(ls-2.76)
152        if (iddist.eq.1) then
153          do ig=1,ngrid
154             topdust(ig)=topdustref         ! constant dust layer top
155          end do
156
157        else if (iddist.eq.2) then          ! "Viking" scenario
158          do ig=1,ngrid
159            topdust(ig)=topdust0(ig)+18.*zlsconst
160          end do
161
162        else if(iddist.eq.3) then         !"MGS" scenario
163          do ig=1,ngrid
164            topdust(ig)=60.+18.*zlsconst
165     &                -(32+18*zlsconst)*sin(lati(ig))**4
166     &                 - 8*zlsconst*(sin(lati(ig)))**5
167          end do
168        endif
169
170
171c       Optical depth in each layer :
172c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
173        if(iddist.ge.1) then
174         expfactor=0.
175          DO l=1,nlayer
176             DO ig=1,ngrid
177             if(pplay(ig,l).gt.700.
178     $                        /(988.**(topdust(ig)/70.))) then
179                zp=(700./pplay(ig,l))**(70./topdust(ig))
180                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
181              else   
182               expfactor=1.e-3
183              endif       
184                aerosol(ig,l,1)= tauref(ig)/700. *
185     s           (pplev(ig,l)-pplev(ig,l+1))
186     &           *expfactor
187c     s           *max( exp(.007*(1.-max(zp,1.))) , 1.E-3 )
188             ENDDO
189          ENDDO
190c     changement dans le calcul de la distribution verticale         
191c     dans le cas des scenarios de poussieres assimiles
192c        if (iaervar.eq.4) THEN  ! TES
193c        call zerophys(ngrid*naerkind,tau)
194c
195c        do l=1,nlayer
196c           do ig=1,ngrid
197c                tau(ig,1)=tau(ig,1)+ aerosol(ig,l,1)
198c           end do
199c        end do
200c        do l=1,nlayer
201c           do ig=1,ngrid
202c               aerosol(ig,l,1)=aerosol(ig,l,1)*tauref(ig)/tau(ig,1)
203c     $     *(pplev(ig,1)/700)
204c           end do
205c        end do
206c        endif
207cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc       
208        else if(iddist.eq.0) then   
209c         old dust vertical distribution function (pollack90)
210          DO l=1,nlayer
211             DO ig=1,ngrid
212                zp=700./pplay(ig,l)
213                aerosol(ig,l,1)= tauref(ig)/700. *
214     s           (pplev(ig,l)-pplev(ig,l+1))
215     s           *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 )
216             ENDDO
217          ENDDO
218        end if
219
220c  ---------------------------------------------------------------------
221c     2) Transported radiatively active dust  (if tracer=T and active=T)
222c ----------------------------------------------------------------------
223      ELSE  IF ((tracer) .and. (active)) THEN
224c     The dust opacity is computed from q
225
226c     a) "doubleq" technique (transport of mass and number mixing ratio)
227c        ~~~~~~~~~~~~~~~~~~~
228       if(doubleq) then
229
230         call zerophys(ngrid*nlayer*naerkind,aerosol)
231
232c        Computing effective radius :
233         do  l=1,nlayer
234           do ig=1, ngrid
235              r0=
236     &        (r3n_q*pq(ig,l,1)/max(pq(ig,l,2),0.01))**(1./3.)
237              r0=min(max(r0,1.e-10),500.e-6)
238              reff= ref_r0 * r0
239cc           If  reff is small, the transported dust mean Qext
240c            is reduced from the reference dust Qext by a factor "coefsize"
241
242             coefsize=min(max(2.52e6*reff-0.043 ,0.)   ,1.)
243
244cc              It is added 1.e-8 to pq to avoid low
245
246                aerosol(ig,l,1)=aerosol(ig,l,1)+ 1.E-8 +
247     &         ( 0.75*Qext(1)*coefsize/(rho_dust*reff))
248     &          * (pq(ig,l,1))*
249     &          (pplev(ig,l)-pplev(ig,l+1))/g
250           end do
251         end do
252         call zerophys(ngrid,tauref)
253
254c     b) Size bin technique (each aerosol can contribute to opacity))
255c        ~~~~~~~~~~~~~~~~~~
256       else
257c        The dust opacity is computed from q
258         call zerophys(ngrid*nlayer*naerkind,aerosol)
259         do iq=1,dustbin
260           do l=1,nlayer
261              do ig=1,ngrid
262cc               qextrhor(iq) is  (3/4)*Qext/(rho*reff)
263cc               It is added 1.e-8 to pq to avoid low
264                 aerosol(ig,l,1)=aerosol(ig,l,1)+
265     &           qextrhor(iq)* (pq(ig,l,iq) + 1.e-8)*
266     &           (pplev(ig,l)-pplev(ig,l+1))/g
267              end do
268           end do
269         end do
270         call zerophys(ngrid,tauref)
271       end if  ! (doubleq)
272      END IF   ! (dust scenario)
273
274c --------------------------------------------------------------------------
275c Column integrated visible optical depth in each point (used for diagnostic)
276c --------------------------------------------------------------------------
277      call zerophys(ngrid*naerkind,tau)
278      do l=1,nlayer
279         do ig=1,ngrid
280               tau(ig,1)=tau(ig,1)+ aerosol(ig,l,1)
281         end do
282      end do
283
284      return
285      end
286
Note: See TracBrowser for help on using the repository browser.