source: trunk/LMDZ.MARS/libf/phymars/thermcell_dqupdown.F90 @ 337

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

MESOSCALE: changes related to g95 support. but the model does not compile apparently. notes taken in NOTES.txt but sorry g95 with LMD_MM_MARS no longer supported at least on 64bits machines. LMDZ.MARS: minor changes related to warning while trying to compile with g95.

  • Property svn:executable set to *
File size: 10.1 KB
Line 
1      subroutine thermcell_dqupdown(ngrid,nlayer,ptimestep,fm0,entr0,  &
2     &    detr0,masse0,q_therm,dq_therm,ztvd,fm_down,ztv,charvar,lmax)
3      implicit none
4
5! #include "iniprint.h"
6!=======================================================================
7!
8!   Calcul du transport verticale dans la couche limite en presence
9!   de "thermiques" explicitement representes
10!   calcul du dq/dt une fois qu'on connait les ascendances
11!   Version modifiee pour prendre les downdrafts a la place de la
12!   subsidence compensatoire
13!=======================================================================
14
15#include "dimensions.h"
16#include "dimphys.h"
17
18! ============================ INPUTS ============================
19
20      INTEGER, INTENT(IN) :: ngrid,nlayer
21      REAL, INTENT(IN) :: ptimestep
22      REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1)
23      REAL, INTENT(IN) :: entr0(ngridmx,nlayermx),detr0(ngridmx,nlayermx)
24      REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
25      REAL, INTENT(IN) :: fm_down(ngridmx,nlayermx+1)
26      REAL, INTENT(IN) :: ztvd(ngridmx,nlayermx)
27      REAL, INTENT(IN) :: ztv(ngridmx,nlayermx)
28      CHARACTER (LEN=20), INTENT(IN) :: charvar
29      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
30      INTEGER, INTENT(IN) :: lmax(ngridmx)
31
32! ============================ OUTPUTS ===========================
33
34      REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)
35
36! ============================ LOCAL =============================
37
38!      REAL detr0(ngridmx,nlayermx)
39      REAL detrd(ngridmx,nlayermx)
40      REAL entrd(ngridmx,nlayermx)     
41      REAL fmd(ngridmx,nlayermx+1)
42      REAL q(ngridmx,nlayermx)
43      REAL qa(ngridmx,nlayermx)
44      REAL qd(ngridmx,nlayermx)
45      INTEGER ig,k
46      LOGICAL active(ngridmx,nlayermx)
47      INTEGER lmax_down(ngridmx),lmin_down(ngridmx)
48      INTEGER ncorec
49
50! =========== Init ==============================================
51
52      entrd(:,:)=0.
53      detrd(:,:)=0.
54      qa(:,:)=q_therm(:,:)
55      q(:,:)=q_therm(:,:)
56!      qd(:,:)=q_therm(:,:)
57!      active(:,:)=.false.
58
59! previous calculation of zdthl_down uses the divergence of fmd
60! so it can be negative without problem. Here we include the sign
61! of fmd in the equations, so it has to be positive
62!
63!      fmd(:,:)=-fm_down(:,:)
64!
65!! ========== Entrainment, Detrainement and Mass =================
66!
67!! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW ===============
68!
69!      do ig=1,ngridmx
70!          if (ztv(ig,nlayermx)-ztvd(ig,nlayermx) .gt. 0.5) then
71!            print*,"downdraft non nul derniere couche !!! (dqupdown)"
72!          endif
73!          detrd(ig,nlayermx)=0.
74!          entrd(ig,nlayermx)=0.
75!      enddo
76!
77!      do k=nlayermx-1,1,-1
78!          do ig=1,ngridmx
79!
80!          if (ztv(ig,k)-ztvd(ig,k) .gt. 0.0001) then
81!          detrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztv(ig,k)-ztvd(ig,k+1)))     &
82!     &    /(ztv(ig,k)-ztvd(ig,k)) - fmd(ig,k))
83!          entrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztvd(ig,k)-ztvd(ig,k+1)))    &
84!     &    /(ztv(ig,k)-ztvd(ig,k)))     
85!
86!          endif
87!        enddo
88!      enddo
89!
90!! ======= We have computed entrainment and detrainment from a prescribed
91!! mass flux and potential temp profile. Due to the way downdraft are parametrized,
92!! this can yield negative entr and detr. We force it to be positive, but in
93!! order to conserve tracers, we need to recompute an adequate mass flux
94!! and modify interface rates, to preserve consistency.
95!      lmax_down(:)=1
96!      lmin_down(:)=1
97!
98!      do k=1,nlayermx
99!        do ig=1,ngridmx
100!        if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
101!!         if (entrd(ig,k).gt.detrd(ig,k)) then
102!        lmax_down(ig)=min(k,lmax(ig))
103!        endif
104!        enddo
105!      enddo
106!      do k=nlayermx,1,-1
107!        do ig=1,ngridmx
108!        if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
109!!         if (detrd(ig,k).gt.entrd(ig,k)) then
110!        lmin_down(ig)=k
111!        endif
112!        enddo
113!      enddo
114!
115!
116!      fmd(:,:)=0.
117!       
118!      do ig=1,ngridmx
119!          if ((lmax_down(ig).gt.1) .and. ((lmax_down(ig)-lmin_down(ig)).gt.1)) then
120!!          fmd(ig,lmax_down(ig))=0.
121!!          entrd(ig,lmax_down(ig))=detrd(ig,lmax_down(ig))
122!!          detrd(ig,lmax_down(ig))=0.
123!!           print*,lmin_down(ig),lmax_down(ig),lmax(ig)
124!
125!            fmd(ig,lmax_down(ig)+1)=0.
126!
127!      do k=lmax_down(ig),lmin_down(ig)+1,-1
128!            fmd(ig,k)=fmd(ig,k+1)+entrd(ig,k)-detrd(ig,k)
129!      enddo
130!           
131!            fmd(ig,lmin_down(ig))=0.
132!            detrd(ig,lmin_down(ig))=fmd(ig,lmin_down(ig)+1)+entrd(ig,lmin_down(ig))
133!
134!          else
135!           entrd(ig,:)=0.
136!           detrd(ig,:)=0.
137!           active(ig,:)=.false.
138!          endif
139!
140!      enddo
141!         ncorec=0
142!         do k=nlayermx,2,-1
143!         do ig=1,ngridmx
144!            if (fmd(ig,k).lt.0.) then
145!!               detrd(ig,k)=max(0.,detrd(ig,k)+fmd(ig,k-1))
146!!               fmd(ig,k-1)=0.
147!!               entrd(ig,k-1)=0.
148!!               detrd(ig,k-1)=0.
149!!               lmin_down(ig)=k-1
150!                fmd(ig,k)=fmd(ig,k+1)
151!                detrd(ig,k)=entrd(ig,k)
152!                ncorec=ncorec+1
153!!                fmd(ig,k)=0.
154!!                detrd(ig,k)=entrd(ig,k)+fmd(ig,k+1)
155!
156!            endif
157!         enddo
158!         enddo
159!
160!         if (ncorec .ne. 0) then
161!         print*, 'corrections for negative downward mass flux :',ncorec
162!         endif
163!         print*, lmin_down(:),lmax_down(:)
164!
165!      do k=2,nlayermx
166!      do ig=1,ngridmx
167!          active(ig,k)=(k.ge.lmin_down(ig)).and.(k.le.lmax_down(ig)) &
168!     & .and.(((fmd(ig,k)+detrd(ig,k))*ptimestep).gt.1.e-6*masse0(ig,k))
169!      enddo
170!      enddo
171!
172!
173!      do ig=1,ngridmx
174!      do k=lmin_down(ig),lmax_down(ig)
175!          if(.not.active(ig,k)) then
176!             active(ig,:)=.false.
177!          endif
178!      enddo
179!      enddo
180!
181!      if(charvar .eq. 'tke') then
182!      active(:,:)=.false.
183!      endif
184!
185!!      do ig=1,ngridmx
186!!         active(ig,lmax_down(ig))=(((fmd(ig,lmax_down(ig))+detrd(ig,lmax_down(ig)))* &
187!!     &       ptimestep).gt.1.e-8*masse0(ig,lmax_down(ig)))
188!!      enddo
189!!
190!!      do ig=1,ngridmx
191!!          if (lmax_down(ig).gt.1) then
192!!            do k=lmax_down(ig)-1,lmin_down(ig),-1
193!!            active(ig,k)=(((fmd(ig,k)+detrd(ig,k))* &
194!!     &       ptimestep).gt.1.e-8*masse0(ig,k)) &
195!!     &     .and. active(ig,k+1)
196!!            enddo
197!!          else
198!!            active(ig,:)=.false.
199!!          endif
200!!      enddo
201!!
202!! ========== qa : q in updraft ==================================
203!
204      do k=2,nlayermx
205         do ig=1,ngridmx
206            if ((fm0(ig,k+1)+detr0(ig,k))*ptimestep.gt.  &
207     &         1.e-5*masse0(ig,k)) then
208         qa(ig,k)=(fm0(ig,k)*qa(ig,k-1)+entr0(ig,k)*q(ig,k))  &
209     &         /(fm0(ig,k+1)+detr0(ig,k))
210
211            if ((qa(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then
212                 if ((q(ig,k) .gt. 0.) .and. ( q(ig,k) .gt. 1.e-15 )) then
213                   !! only print if it is the thermal scheme which makes qa<0
214                   print*,'qa<0 created by thermals!!!',charvar,ig,k, &
215     &                     fm0(ig,k),qa(ig,k-1),     &
216     &                     entr0(ig,k),q(ig,k),fm0(ig,k+1),detr0(ig,k)
217                   print*,'---------> Cancelling qa'
218                 endif
219                 qa(ig,k)=q(ig,k)  !! no action of thermals in this case.
220            endif
221            else
222               qa(ig,k)=q(ig,k)
223            endif
224            enddo
225      enddo
226
227! ========== qd : q in downdraft =================================
228!
229!      do k=nlayermx-1,1,-1
230!         do ig=1,ngridmx
231!            if (active(ig,k)) then
232!         qd(ig,k)=(fmd(ig,k+1)*qd(ig,k+1)+entrd(ig,k)*q(ig,k))  &
233!     &         /(fmd(ig,k)+detrd(ig,k))
234!
235!            if ((qd(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then
236!     print*,'qd<0!!!',charvar,ig,k,fmd(ig,k),qd(ig,k),entrd(ig,k),q(ig,k),fmd(ig,k+1),detrd(ig,k),lmin_down(ig),lmax_down(ig)
237!     print*, '---------> cancelling qd, no downdraft for this gridpoint'
238!                     qd(ig,k)=q(ig,k)
239!                     active(ig,:)=.false.
240!            endif
241!            else
242!               qd(ig,k)=q(ig,k)
243!            endif
244!!          print*,'active,k,entr,detr,q,qd (down) :',active(ig,k),k,entrd(ig,k),detrd(ig,k),q(ig,k),qd(ig,k)
245!         enddo
246!      enddo
247!
248!
249! ====== dq ======================================================
250
251      do ig=1,ngridmx
252!         if(active(ig,1)) then
253!
254!         dq_therm(ig,1)=(detr0(ig,1)*qa(ig,1)+detrd(ig,1)*qd(ig,1) &
255!      &               +fm0(ig,2)*q(ig,2)   &
256!      &               -entr0(ig,1)*q(ig,1)-entrd(ig,1)*q(ig,1)   &
257!      &               -fmd(ig,2)*q(ig,1)) &
258!      &               *ptimestep/masse0(ig,1)
259!
260!         else
261         dq_therm(ig,1)=(detr0(ig,1)*qa(ig,1)+fm0(ig,2)*q(ig,2) &
262      &               -entr0(ig,1)*q(ig,1)) &
263      &               *ptimestep/masse0(ig,1)
264
265!         endif
266       enddo
267     
268       do k=2,nlayermx-1
269         do ig=1, ngridmx
270
271!         if(active(ig,k)) then
272!
273!         dq_therm(ig,k)=(detr0(ig,k)*qa(ig,k)+detrd(ig,k)*qd(ig,k) &
274!      &               +fm0(ig,k+1)*q(ig,k+1)+fmd(ig,k)*q(ig,k-1)   &
275!      &               -entr0(ig,k)*q(ig,k)-entrd(ig,k)*q(ig,k)   &
276!      &               -fm0(ig,k)*q(ig,k)-fmd(ig,k+1)*q(ig,k))      &
277!     &               *ptimestep/masse0(ig,k)
278
279
280!         else
281         dq_therm(ig,k)=(detr0(ig,k)*qa(ig,k)+fm0(ig,k+1)*q(ig,k+1) &
282      &               -entr0(ig,k)*q(ig,k)-fm0(ig,k)*q(ig,k))  &
283      &               *ptimestep/masse0(ig,k)
284
285
286!         endif
287
288         enddo
289      enddo
290
291         do ig=1, ngridmx
292
293!         if(active(ig,nlayermx)) then
294!
295!         dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx)+detrd(ig,nlayermx)*qd(ig,nlayermx) &
296!      &               +fmd(ig,nlayermx)*q(ig,nlayermx-1)   &
297!      &         -entr0(ig,nlayermx)*q(ig,nlayermx)-entrd(ig,nlayermx)*q(ig,nlayermx)   &
298!      &               -fm0(ig,nlayermx)*q(ig,nlayermx)) &
299!      &               *ptimestep/masse0(ig,nlayermx)
300
301!         else
302         dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx) &
303      &             -entr0(ig,nlayermx)*q(ig,nlayermx)-fm0(ig,nlayermx)*q(ig,nlayermx)) &
304      &               *ptimestep/masse0(ig,nlayermx)
305!         endif
306         
307         enddo
308      return
309      end
Note: See TracBrowser for help on using the repository browser.