source: trunk/LMDZ.MARS/libf/phymars/dustdevil.F @ 171

Last change on this file since 171 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM
File size: 7.0 KB
Line 
1      SUBROUTINE dustdevil(ngrid,nlay,nq, pplev,pu,pv,pt, ptsurf,pq2,
2     &                pdqdev,pdqs_dev)
3      IMPLICIT NONE
4
5c=======================================================================
6c
7c
8c  VERSION SPECIAL TRACEURS :
9c  Parameterization of dust devil activities
10c  to estimate dust lifting
11c  The dust devil activity is estimated based on
12c  Renno et al. 1998 (JAS 55, 3244-3252) 
13c
14c  It is proportional to (1-b)*Fs
15c
16c  With b= [ps**(rcp+1) - ptop**(rcp+1)] / [(ps-ptop)*(rcp+1)* ps**rcp]
17c  with ptop pressure of the top of the boundary layer
18c       rcp = R/cp
19c
20c  and Fs the surface sensible heat flux = Cd*|U|*(T(1) -Tsurf)
21c       
22c=======================================================================
23
24c-----------------------------------------------------------------------
25c   declarations:
26c   -------------
27
28#include "dimensions.h"
29#include "dimphys.h"
30#include "comcstfi.h"     
31c#include "comconst.h"        ! TEMPORAIRE AVEC ANLDEVIL !!!!
32#include "planete.h"
33#include "comgeomfi.h"
34#include "tracer.h"
35c   arguments:
36c   ----------
37
38      INTEGER ngrid,nlay
39      REAL pplev(ngrid,nlay+1)
40      REAL pt(ngrid,nlay)
41      REAL pu(ngrid,nlay)
42      REAL pv(ngrid,nlay)
43      REAL pq2(ngrid,nlay+1)
44      REAL ptsurf(ngrid)
45
46c    Traceurs :
47      integer nq
48      real pdqdev(ngrid,nlay,nq)
49      real pdqs_dev(ngrid,nq)
50     
51c   local:
52c   ------
53
54      INTEGER ig,l,iq
55      real Cd, z1
56      save Cd
57
58      LOGICAL firstcall
59      SAVE firstcall
60
61
62      REAL devila(ngridmx)
63      integer ltop(ngridmx)
64      real b,rho,Fs,wind
65
66
67
68      REAL  q2top , seuil
69      SAVE  q2top , seuil
70      DATA q2top/.5/ ! value of q2 at the top of PBL
71      DATA seuil/.3/ ! value of minimum dust devil activity for dust lifting
72
73
74      DATA firstcall/.true./
75
76c   TEMPORAIRE AVEC ANLDEVIL : *************
77c        real b_diag(ngridmx)
78c       real localtime(ngridmx)
79c       common/temporaire/localtime
80c      real ztop(ngridmx),magwind(ngridmx),t1(ngridmx)
81c      real rcp ,cpp
82c      rcp = kappa
83c      cpp = r/rcp
84c   **********************************
85       
86
87c-----------------------------------------------------------------------
88c    initialisation
89c    --------------
90
91      IF (firstcall) THEN
92
93        write(*,*) 'In dustdevil :'
94        write(*,*) '    q2top= ',q2top,'     seuil= ', seuil
95c       un petit test de coherence:
96         IF(ngrid.NE.ngridmx) THEN
97            PRINT*,'STOP dans coefdifv'
98            PRINT*,'probleme de dimensions :'
99            PRINT*,'ngrid  =',ngrid
100            PRINT*,'ngridmx  =',ngridmx
101            STOP
102         ENDIF
103
104c A rough estimation of the horizontal drag coefficient Cd
105c (the scale heigh is taken to be 13 km since we are typically
106c dealing with daytime temperature around 250K.
107c
108         z1= -0.5*13.e3*log(pplev(1,2)/pplev(1,1))
109         Cd = (0.4/log(z1/z0))**2
110
111         firstcall=.false.
112
113c        Temporaire
114c        open(77,file='devil')
115     
116      ENDIF
117
118c-----------------------------------------------------------------------
119c Initialisation
120      do iq=1,nq
121       do l=1,nlay
122           do ig=1,ngrid
123             pdqdev(ig,l,iq)= 0
124           end do
125       end do
126      end do
127
128
129c-----------------------------------------------------------------------
130c      Determining the top of the boundary layer
131c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132      do ig=1,ngrid
133         do  l=2,nlay-1
134            if (pq2(ig,l).lt.q2top)then
135              ltop(ig)=l
136              goto 99
137            end if
138         enddo
139 99      continue
140
141c        ***************************************
142cc        if (ptsurf(ig).gt.255)then
143c         write(*,*) 'tsurf, ztop (km): ', ig,
144c     &   ptsurf(ig), -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
145c         if ((ptsurf(ig).gt.50.).and.(
146c     &      (-12.*log(pplev(ig,ltop(ig))/pplev(ig,1))).gt.60.))then
147c            do l=1,nlay
148c             write(*,*) l,
149c     &       -12.*log(pplev(ig,l)/pplev(ig,1)),pq2(ig,l)
150c            end do
151c         end if
152cc        end if
153c        ***************************************
154     
155      enddo
156
157c        ***************************************
158c        do ig=100,148
159c           write(*,*)'time,tsurf,ztop', localtime(ig),ptsurf(ig),
160c    &      -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
161c        end do
162c        ***************************************
163
164
165c   Calculation : dust devil intensity
166c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
167      do ig=1,ngrid
168
169c --------------------------------------------------
170c 1) Version 1 : sensible heat flux using actual wind :
171c        Wind magnitude:
172c        wind = sqrt(pu(ig,1)**2+pv(ig,1)**2)
173c
174c --------------------------------------------------
175c 2) Version 2 : sensible heat flux using  wind = 15 m/s
176         wind = 15.
177c ----------------------------------------------------
178c        Density :
179         rho=pplev(ig,1)/(R*pt(ig,1))
180
181c        Sensible heat flux (W.m-2) (>0 if up)
182         Fs= rho*cpp*Cd * wind
183     &       * (ptsurf(ig) -pt(ig,1))
184         b= (pplev(ig,1)**(rcp+1) - pplev(ig,ltop(ig))**(rcp+1)) /
185     &    ( (pplev(ig,1)-pplev(ig,ltop(ig)))*(rcp+1)*pplev(ig,1)**rcp)
186
187c        b_diag(ig) = b     ! TEMPORAIRE (pour diagnostique)
188
189c   Energy flux available to drive dust devil (W.m-2) : (1-b)*Fs
190c   Dust devil activity :
191         devila(ig)= max( 0. , (1-b)*Fs - seuil )
192      enddo
193c   
194c     lifted dust (kg m-2 s-1)  (<0 when lifting)
195c     ~~~~~~~~~~ 
196      do iq=1,nq
197         do ig=1,ngrid
198           pdqs_dev(ig,iq)= - alpha_devil(iq) * devila(ig)
199         end do
200      end do
201
202c     Injection of dust in the atmosphere (up to the top of pbl)
203c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
204      do iq=1,nq
205       do ig=1,ngrid
206         if (devila(ig).ne.0.) then
207           do l=1,ltop(ig)
208             pdqdev(ig,l,iq)=-pdqs_dev(ig,iq)*g/
209     &        (pplev(ig,1)-pplev(ig,ltop(ig)))
210           end do
211         end if
212       end do
213      end do
214
215c *********************************************************     
216c     TEMPORAIRE AVEC ANLDEVIL:
217c     IF (ngrid.gt.1) THEN
218c      do ig=2,ngridmx-1
219c       write(77,88) lati(ig)*180./pi,localtime(ig),
220c    &        -12.*log(pplev(ig,ltop(ig))/pplev(ig,1)),
221c    &   devil(ig),min(sqrt(pu(ig,1)**2+pv(ig,1)**2),40.),
222c    &   ptsurf(ig)-pt(ig,1),ptsurf(ig),pplev(ig,1)
223c      end do   
224c88    format (f7.3,1x,f7.3,1x,f6.3,1x,f6.4,1x,f7.4,1x,
225c    &        f7.3,1x,f7.3,1x,f9.3)
226c      do ig=1,ngridmx
227c       ztop(ig) = -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
228c       magwind(ig) = sqrt(pu(ig,1)**2+pv(ig,1)**2)
229c       t1(ig) =max(ptsurf(ig)- pt(ig,1),0.)
230c      end do
231
232c       call WRITEDIAGFI(ngridmx,'dqs_dev','dqs devil',
233c    &               'kg.m-2.s-1',2,pdqs_dev)
234c       call WRITEDIAGFI(ngridmx,'wind','wind',
235c    &               'm.s-1',2,magwind)
236c       call WRITEDIAGFI(ngridmx,'ztop','top pbl',
237c    &               'km',2,ztop)
238c       call WRITEDIAGFI(ngridmx,'tsurf','tsurf',
239c    &               'K',2,ptsurf)
240c       call WRITEDIAGFI(ngridmx,'T1','T(1)',
241c    &               'K',2,t1)
242c       call WRITEDIAGFI(ngridmx,'b','b',
243c    &               ' ',2,b_diag)
244c     END If
245c *********************************************************     
246         
247      RETURN
248      END
249
250
Note: See TracBrowser for help on using the repository browser.