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

Last change on this file since 2616 was 2608, checked in by romain.vande, 3 years ago

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