source: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/thermcell_dq.F90 @ 2279

Last change on this file since 2279 was 983, checked in by Laurent Fairhead, 16 years ago

Probleme d'initialisation YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.0 KB
Line 
1      subroutine thermcell_dq(ngrid,nlay,ptimestep,fm,entr,  &
2     &           masse,q,dq,qa,lev_out)
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!
12!=======================================================================
13
14      integer ngrid,nlay
15
16      real ptimestep
17      real masse(ngrid,nlay),fm(ngrid,nlay+1)
18      real entr(ngrid,nlay)
19      real q(ngrid,nlay)
20      real dq(ngrid,nlay)
21      integer lev_out                           ! niveau pour les print
22
23      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
24
25      real zzm
26
27      integer ig,k
28      real cfl
29
30      real qold(ngrid,nlay)
31      real ztimestep
32      integer niter,iter
33
34
35
36! Calcul du critere CFL pour l'advection dans la subsidence
37      cfl = 0.
38      do k=1,nlay
39         do ig=1,ngrid
40            zzm=masse(ig,k)/ptimestep
41            cfl=max(cfl,fm(ig,k)/zzm)
42            if (entr(ig,k).gt.zzm) then
43               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
44               stop
45            endif
46         enddo
47      enddo
48
49!IM 090508     print*,'CFL CFL CFL CFL ',cfl
50
51#undef CFL
52#ifdef CFL
53! On subdivise le calcul en niter pas de temps.
54      niter=int(cfl)+1
55#else
56      niter=1
57#endif
58
59      ztimestep=ptimestep/niter
60      qold=q
61
62
63do iter=1,niter
64      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
65
66!   calcul du detrainement
67      do k=1,nlay
68         do ig=1,ngrid
69            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
70!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
71!test
72            if (detr(ig,k).lt.0.) then
73               entr(ig,k)=entr(ig,k)-detr(ig,k)
74               detr(ig,k)=0.
75!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
76!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
77            endif
78            if (fm(ig,k+1).lt.0.) then
79!               print*,'fm2<0!!!'
80            endif
81            if (entr(ig,k).lt.0.) then
82!               print*,'entr2<0!!!'
83            endif
84         enddo
85      enddo
86
87!   calcul de la valeur dans les ascendances
88      do ig=1,ngrid
89         qa(ig,1)=q(ig,1)
90      enddo
91
92      do k=2,nlay
93         do ig=1,ngrid
94            if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
95     &         1.e-5*masse(ig,k)) then
96         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
97     &         /(fm(ig,k+1)+detr(ig,k))
98            else
99               qa(ig,k)=q(ig,k)
100            endif
101            if (qa(ig,k).lt.0.) then
102!               print*,'qa<0!!!'
103            endif
104            if (q(ig,k).lt.0.) then
105!               print*,'q<0!!!'
106            endif
107         enddo
108      enddo
109
110! Calcul du flux subsident
111
112      do k=2,nlay
113         do ig=1,ngrid
114#undef centre
115#ifdef centre
116             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
117#else
118
119#define plusqueun
120#ifdef plusqueun
121! Schema avec advection sur plus qu'une maille.
122            zzm=masse(ig,k)/ztimestep
123            if (fm(ig,k)>zzm) then
124               wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
125            else
126               wqd(ig,k)=fm(ig,k)*q(ig,k)
127            endif
128#else
129            wqd(ig,k)=fm(ig,k)*q(ig,k)
130#endif
131#endif
132
133            if (wqd(ig,k).lt.0.) then
134!               print*,'wqd<0!!!'
135            endif
136         enddo
137      enddo
138      do ig=1,ngrid
139         wqd(ig,1)=0.
140         wqd(ig,nlay+1)=0.
141      enddo
142     
143
144! Calcul des tendances
145      do k=1,nlay
146         do ig=1,ngrid
147            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
148     &               -wqd(ig,k)+wqd(ig,k+1))  &
149     &               *ztimestep/masse(ig,k)
150!            if (dq(ig,k).lt.0.) then
151!               print*,'dq<0!!!'
152!            endif
153         enddo
154      enddo
155
156
157enddo
158
159
160! Calcul des tendances
161      do k=1,nlay
162         do ig=1,ngrid
163            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
164            q(ig,k)=qold(ig,k)
165         enddo
166      enddo
167
168      return
169      end
Note: See TracBrowser for help on using the repository browser.