source: LMDZ4/trunk/libf/phylmd/thermcell_dv2.F90 @ 1934

Last change on this file since 1934 was 1408, checked in by Laurent Fairhead, 14 years ago

Ooops, sorry

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.7 KB
RevLine 
[878]1      subroutine thermcell_dv2(ngrid,nlay,ptimestep,fm,entr,masse  &
2     &    ,fraca,larga  &
3     &    ,u,v,du,dv,ua,va,lev_out)
4      implicit none
5
[1408]6#include "iniprint.h"
[878]7!=======================================================================
8!
9!   Calcul du transport verticale dans la couche limite en presence
10!   de "thermiques" explicitement representes
11!   calcul du dq/dt une fois qu'on connait les ascendances
12!
[1403]13! Vectorisation, FH : 2010/03/08
14!
[878]15!=======================================================================
16
17
18      integer ngrid,nlay
19
20      real ptimestep
21      real masse(ngrid,nlay),fm(ngrid,nlay+1)
22      real fraca(ngrid,nlay+1)
23      real larga(ngrid)
24      real entr(ngrid,nlay)
25      real u(ngrid,nlay)
26      real ua(ngrid,nlay)
27      real du(ngrid,nlay)
28      real v(ngrid,nlay)
29      real va(ngrid,nlay)
30      real dv(ngrid,nlay)
31      integer lev_out                           ! niveau pour les print
32
33      real qa(ngrid,nlay),detr(ngrid,nlay),zf,zf2
34      real wvd(ngrid,nlay+1),wud(ngrid,nlay+1)
[1403]35      real gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1)
[878]36      real ue(ngrid,nlay),ve(ngrid,nlay)
[1403]37      LOGICAL ltherm(ngrid,nlay)
38      real dua(ngrid,nlay),dva(ngrid,nlay)
[878]39      integer iter
40
[1403]41      integer ig,k,nlarga0
[878]42
[1403]43!-------------------------------------------------------------------------
44
[878]45!   calcul du detrainement
[1403]46!---------------------------
[878]47
[1407]48!      print*,'THERMCELL DV2 OPTIMISE 3'
[1403]49
50      nlarga0=0.
51
[878]52      do k=1,nlay
53         do ig=1,ngrid
54            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
55         enddo
56      enddo
57
58!   calcul de la valeur dans les ascendances
59      do ig=1,ngrid
60         ua(ig,1)=u(ig,1)
61         va(ig,1)=v(ig,1)
62         ue(ig,1)=u(ig,1)
63         ve(ig,1)=v(ig,1)
64      enddo
65
[1146]66      IF(prt_level>9)WRITE(lunout,*)                                    &
67     &      'WARNING on initialise gamma(1:ngrid,1)=0.'
[972]68      gamma(1:ngrid,1)=0.
[878]69      do k=2,nlay
70         do ig=1,ngrid
[1403]71            ltherm(ig,k)=(fm(ig,k+1)+detr(ig,k))*ptimestep > 1.e-5*masse(ig,k)
72            if(ltherm(ig,k).and.larga(ig)>0.) then
73               gamma0(ig,k)=masse(ig,k)  &
[878]74     &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
75     &         *0.5/larga(ig)  &
76     &         *1.
[1403]77            else
78               gamma0(ig,k)=0.
79            endif
80            if (ltherm(ig,k).and.larga(ig)<=0.) nlarga0=nlarga0+1
81         enddo
82      enddo
83
84      gamma(:,:)=0.
85
86      do k=2,nlay
87
88         do ig=1,ngrid
89            if (ltherm(ig,k)) then
90               dua(ig,k)=ua(ig,k-1)-u(ig,k-1)
91               dva(ig,k)=va(ig,k-1)-v(ig,k-1)
92            else
93               ua(ig,k)=u(ig,k)
94               va(ig,k)=v(ig,k)
95               ue(ig,k)=u(ig,k)
96               ve(ig,k)=v(ig,k)
97            endif
98         enddo
99
100
101! Debut des iterations
102!----------------------
103do iter=1,5
104         do ig=1,ngrid
105! Pour memoire : calcul prenant en compte la fraction reelle
106!              zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
107!              zf2=1./(1.-zf)
108! Calcul avec fraction infiniement petite
[878]109               zf=0.
[1403]110               zf2=1.
111
112!  la première fois on multiplie le coefficient de freinage
113!  par le module du vent dans la couche en dessous.
114!  Mais pourquoi donc ???
115               if (ltherm(ig,k)) then
[878]116!   On choisit une relaxation lineaire.
[1403]117!                 gamma(ig,k)=gamma0(ig,k)
[878]118!   On choisit une relaxation quadratique.
[1403]119                  gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2)
[878]120                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)  &
121     &               +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))  &
122     &               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2  &
123     &                 +gamma(ig,k))
124                  va(ig,k)=(fm(ig,k)*va(ig,k-1)  &
125     &               +(zf2*entr(ig,k)+gamma(ig,k))*v(ig,k))  &
126     &               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2  &
127     &                 +gamma(ig,k))
[1403]128!                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua(ig,k),dva(ig,k)
129                  dua(ig,k)=ua(ig,k)-u(ig,k)
130                  dva(ig,k)=va(ig,k)-v(ig,k)
[878]131                  ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
132                  ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
[1403]133               endif
[878]134         enddo
[1403]135! Fin des iterations
136!--------------------
137enddo
[878]138
[1403]139      enddo ! k=2,nlay
140
141
142! Calcul du flux vertical de moment dans l'environnement.
143!---------------------------------------------------------
[878]144      do k=2,nlay
145         do ig=1,ngrid
146            wud(ig,k)=fm(ig,k)*ue(ig,k)
147            wvd(ig,k)=fm(ig,k)*ve(ig,k)
148         enddo
149      enddo
150      do ig=1,ngrid
151         wud(ig,1)=0.
152         wud(ig,nlay+1)=0.
153         wvd(ig,1)=0.
154         wvd(ig,nlay+1)=0.
155      enddo
156
[1403]157! calcul des tendances.
158!-----------------------
[878]159      do k=1,nlay
160         do ig=1,ngrid
161            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)  &
162     &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
163     &               -wud(ig,k)+wud(ig,k+1))  &
164     &               /masse(ig,k)
165            dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)  &
166     &               -(entr(ig,k)+gamma(ig,k))*ve(ig,k)  &
167     &               -wvd(ig,k)+wvd(ig,k+1))  &
168     &               /masse(ig,k)
169         enddo
170      enddo
171
[1403]172
173! Sorties eventuelles.
174!----------------------
175
176   if(prt_level.GE.10) then
177      do k=1,nlay
178         do ig=1,ngrid
179           print*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), &
180     &   entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), &
181     &   masse(ig,k)
182         enddo
183      enddo
184   endif
185!
186     if (nlarga0>0) then
187          print*,'WARNING !!!!!! DANS THERMCELL_DV2 '
188          print*,nlarga0,' points pour lesquels laraga=0. dans un thermique'
189          print*,'Il faudrait decortiquer ces points'
190     endif
191
[878]192      return
193      end
Note: See TracBrowser for help on using the repository browser.