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

Last change on this file since 2266 was 1779, checked in by aslmd, 7 years ago

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