source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_dq.F90 @ 3722

Last change on this file since 3722 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

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