My Project
 All Classes Files Functions Variables Macros
thermcell_dv2.F90
Go to the documentation of this file.
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 
6 #include "iniprint.h"
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 !
13 ! Vectorisation, FH : 2010/03/08
14 !
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)
35  real gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1)
36  real ue(ngrid,nlay),ve(ngrid,nlay)
37  LOGICAL ltherm(ngrid,nlay)
38  real dua(ngrid,nlay),dva(ngrid,nlay)
39  integer iter
40 
41  integer ig,k,nlarga0
42 
43 !-------------------------------------------------------------------------
44 
45 ! calcul du detrainement
46 !---------------------------
47 
48 ! print*,'THERMCELL DV2 OPTIMISE 3'
49 
50  nlarga0=0.
51 
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 
66  IF(prt_level>9)WRITE(lunout,*) &
67  & 'WARNING on initialise gamma(1:ngrid,1)=0.'
68  gamma(1:ngrid,1)=0.
69  do k=2,nlay
70  do ig=1,ngrid
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) &
74  & *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) ) &
75  & *0.5/larga(ig) &
76  & *1.
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 !----------------------
103 do 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
109  zf=0.
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
116 ! On choisit une relaxation lineaire.
117 ! gamma(ig,k)=gamma0(ig,k)
118 ! On choisit une relaxation quadratique.
119  gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2)
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))
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)
131  ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
132  ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
133  endif
134  enddo
135 ! Fin des iterations
136 !--------------------
137 enddo
138 
139  enddo ! k=2,nlay
140 
141 
142 ! Calcul du flux vertical de moment dans l'environnement.
143 !---------------------------------------------------------
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 
157 ! calcul des tendances.
158 !-----------------------
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 
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 
192  return
193  end