source: LMDZ6/branches/Portage_acc/libf/phylmdiso/thermcell_down.F90 @ 4500

Last change on this file since 4500 was 4447, checked in by Laurent Fairhead, 17 months ago

Added some routines from the trunk that were previously links and that svn did not want to commit on the previous commit (with an
"Node filename has unexpectedly changed kind" error)
File size: 10.4 KB
Line 
1SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,trac,dtrac)
2
3USE thermcell_ini_mod, ONLY: iflag_thermals_down
4
5
6!-----------------------------------------------------------------
7! thermcell_updown_dq: computes the tendency of tracers associated
8! with the presence of convective up/down drafts
9! This routine that has been collectively written during the
10! "ateliers downdrafts" in 2022/2023
11! Maelle, Frédéric, Catherine, Fleur, Florent, Etienne
12!------------------------------------------------------------------
13
14
15   IMPLICIT NONE
16
17! declarations
18!==============================================================
19
20! input/output
21
22   integer,intent(in)  :: ngrid ! number of horizontal grid points
23   integer, intent(in) :: nlay  ! number of vertical layers
24   real,intent(in) :: ptimestep ! time step of the physics [s]
25   real,intent(in), dimension(ngrid,nlay) :: eup ! entrainment to updrafts * dz [same unit as flux]
26   real,intent(in), dimension(ngrid,nlay) :: dup ! detrainment from updrafts * dz [same unit as flux]
27   real,intent(in), dimension(ngrid,nlay) :: edn ! entrainment to downdrafts * dz [same unit as flux]
28   real,intent(in), dimension(ngrid,nlay) :: ddn ! detrainment from downdrafts * dz [same unit as flux]
29   real,intent(in), dimension(ngrid,nlay) :: masse ! mass of layers = rho dz
30   real,intent(in), dimension(ngrid,nlay) :: trac ! tracer
31   integer, intent(in), dimension(ngrid) :: lmax ! max level index at which downdraft are present
32   real,intent(out),dimension(ngrid,nlay) ::dtrac ! tendance du traceur
33
34   
35! Local
36
37   real, dimension(ngrid,nlay+1) :: fup,fdn,fc,fthu,fthd,fthe,fthtot
38   real, dimension(ngrid,nlay) :: tracu,tracd,traci,tracold
39   real :: www, mstar_inv
40   integer ig,ilay
41   real, dimension(ngrid,nlay):: s1,s2,num !coefficients pour la resolution implicite
42   integer :: iflag_impl=1 ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement
43   
44   fdn(:,:)=0.
45   fup(:,:)=0.
46   fc(:,:)=0.
47   fthu(:,:)=0.
48   fthd(:,:)=0.
49   fthe(:,:)=0.
50   fthtot(:,:)=0.
51   tracd(:,:)=0.
52   tracu(:,:)=0.
53   traci(:,:)=trac(:,:)
54   tracold(:,:)=trac(:,:)
55   s1(:,:)=0.
56   s2(:,:)=0.
57   num(:,:)=1.
58
59   if ( iflag_thermals_down < 10 ) then
60        stop 'thermcell_down_dq = 0 or >= 10'
61   else
62        iflag_impl=iflag_thermals_down-10
63   endif
64     
65
66   ! lmax : indice tel que fu(kmax+1)=0
67   ! Dans ce cas, pas besoin d'initialiser tracd(lmax) ( =trac(lmax) )
68   ! Boucle pour le downdraft
69   do ilay=nlay,1,-1
70      do ig=1,ngrid
71         !if ( lmax(ig) > nlay - 2 ) stop "les thermiques montent trop haut"
72         if (ilay.le.lmax(ig) .and. lmax(ig)>1 ) then
73            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
74            if ( fdn(ig,ilay)+ddn(ig,ilay) > 0. ) then
75               www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
76            else
77               www=0.
78            endif
79            tracd(ig,ilay)=www*tracd(ig,ilay+1) + (1.-www)*trac(ig,ilay)
80         endif
81      enddo
82   enddo !Fin boucle sur l'updraft
83   fdn(:,1)=0.
84
85   !Boucle pour l'updraft
86   do ilay=1,nlay,1
87      do ig=1,ngrid
88         if (ilay.lt.lmax(ig) .and. lmax(ig)>1) then
89            fup(ig,ilay+1)=fup(ig,ilay)+eup(ig,ilay)-dup(ig,ilay)
90            if (fup(ig,ilay+1)+dup(ig,ilay) > 0.) then
91               www=fup(ig,ilay)/(fup(ig,ilay+1)+dup(ig,ilay))
92            else
93               www=0.
94            endif
95            if (ilay == 1 ) then
96               tracu(ig,ilay)=trac(ig,ilay)
97            else
98               tracu(ig,ilay)=www*tracu(ig,ilay-1)+(1.-www)*trac(ig,ilay)
99            endif
100         endif
101      enddo
102      enddo !fin boucle sur le downdraft
103
104   ! Calcul des flux des traceurs dans les updraft et les downdrfat
105   ! et du flux de masse compensateur
106   ! en ilay=1 et nlay+1, fthu=0 et fthd=0
107   fthu(:,1)=0.
108   fthu(:,nlay+1)=0.
109   fthd(:,1)=0.
110   fthd(:,nlay+1)=0.
111   fc(:,1)=0.
112   fc(:,nlay+1)=0.
113   do ilay=2,nlay,1 !boucle sur les interfaces
114     do ig=1,ngrid
115       fthu(ig,ilay)=fup(ig,ilay)*tracu(ig,ilay-1)
116       fthd(ig,ilay)=-fdn(ig,ilay)*tracd(ig,ilay)
117       fc(ig,ilay)=fup(ig,ilay)-fdn(ig,ilay)
118     enddo
119   enddo
120   
121
122   !Boucle pour calculer le flux du traceur flux updraft, flux downdraft, flux compensatoire
123   !Methode explicite :
124   if(iflag_impl==0) then
125     do ilay=2,nlay,1
126       do ig=1,ngrid
127         !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher trac au dessus!!!!!
128         !!!! tentative de prise en compte d'un flux compensatoire montant  !!!!
129         if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
130            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
131            stop         
132            !fthe(ig,ilay)=(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1)
133         else
134            fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
135         endif
136         !! si on voulait le prendre en compte on
137         !fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1)
138         fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
139       enddo
140     enddo
141     !Boucle pour calculer trac
142     do ilay=1,nlay
143       do ig=1,ngrid
144         dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))/masse(ig,ilay)
145!         trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
146       enddo
147     enddo !fin du calculer de la tendance du traceur avec la methode explicite
148
149   !!! Reecriture du schéma explicite avec les notations du schéma implicite
150   else if(iflag_impl==-1) then
151     write(*,*) 'nouveau schéma explicite !!!'
152     !!! Calcul de s1
153     do ilay=1,nlay
154       do ig=1,ngrid
155         s1(ig,ilay)=fthu(ig,ilay)-fthu(ig,ilay+1)+fthd(ig,ilay)-fthd(ig,ilay+1)
156         s2(ig,ilay)=s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1)
157       enddo
158     enddo
159
160     do ilay=2,nlay,1
161       do ig=1,ngrid
162         if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
163            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
164            stop         
165         else
166            fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
167         endif
168         fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
169       enddo
170     enddo
171     !Boucle pour calculer trac
172     do ilay=1,nlay
173       do ig=1,ngrid
174         ! dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))/masse(ig,ilay)
175         dtrac(ig,ilay)=(s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1))/masse(ig,ilay)
176!         trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
177!         trac(ig,ilay)=trac(ig,ilay) + (s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1))*(ptimestep/masse(ig,ilay))
178       enddo
179     enddo !fin du calculer de la tendance du traceur avec la methode explicite
180
181   else if (iflag_impl==1) then
182     do ilay=1,nlay
183       do ig=1,ngrid
184         s1(ig,ilay)=fthu(ig,ilay)-fthu(ig,ilay+1)+fthd(ig,ilay)-fthd(ig,ilay+1)
185       enddo
186     enddo
187     
188     !Boucle pour calculer traci = trac((t+dt)
189     do ilay=nlay-1,1,-1
190       do ig=1,ngrid
191         if((fup(ig,ilay)-fdn(ig,ilay)) .lt. 0) then
192            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq dans le cas d une resolution implicite, ilay : ', ilay
193            stop 
194         else
195           mstar_inv=ptimestep/masse(ig,ilay)
196           traci(ig,ilay)=((traci(ig,ilay+1)*fc(ig,ilay+1)+s1(ig,ilay))*mstar_inv+tracold(ig,ilay))/(1.+fc(ig,ilay)*mstar_inv)
197         endif
198       enddo
199     enddo
200     do ilay=1,nlay
201       do ig=1,ngrid
202         dtrac(ig,ilay)=(traci(ig,ilay)-tracold(ig,ilay))/ptimestep
203       enddo
204     enddo
205
206   else
207      write(*,*) 'valeur de iflag_impl non prevue'
208      stop
209
210   endif
211
212 RETURN
213   END
214
215!=========================================================================
216
217   SUBROUTINE thermcell_down(ngrid,nlay,po,pt,pu,pv,pplay,pplev,  &
218     &           lmax,fup,eup,dup,theta)
219
220!--------------------------------------------------------------
221!thermcell_down: calcul des propri??t??s du panache descendant.
222!--------------------------------------------------------------
223
224
225   USE thermcell_ini_mod, ONLY : prt_level,RLvCp,RKAPPA,RETV,fact_thermals_down
226   IMPLICIT NONE
227
228! arguments
229
230   integer,intent(in) :: ngrid,nlay
231   real,intent(in), dimension(ngrid,nlay) :: po,pt,pu,pv,pplay,eup,dup
232   real,intent(in), dimension(ngrid,nlay) :: theta
233   real,intent(in), dimension(ngrid,nlay+1) :: pplev,fup
234   integer, intent(in), dimension(ngrid) :: lmax
235
236
237   
238! Local
239
240   real, dimension(ngrid,nlay) :: edn,ddn,thetad
241   real, dimension(ngrid,nlay+1) :: fdn
242
243   integer ig,ilay
244   real dqsat_dT
245   logical mask(ngrid,nlay)
246
247   edn(:,:)=0.
248   ddn(:,:)=0.
249   fdn(:,:)=0.
250   thetad(:,:)=0.
251
252   ! lmax : indice tel que fu(kmax+1)=0
253   
254   ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) )
255
256! FH MODIFS APRES REUNIONS POUR COMMISSIONS
257! quelques erreurs de declaration
258! probleme si lmax=1 ce qui a l'air d'??tre le cas en d??but de simu. Devrait ??tre 0 ?
259! Remarques :
260! on pourrait ??crire la formule de thetad
261!    www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
262!    thetad(ig,ilay)= www * thetad(ig,ilay+1) + (1.-www) * theta(ig,ilay)
263! Elle a l'avantage de bien montr?? la conservation, l'id??e fondamentale dans le
264!   transport qu'on ne fait que sommer des "sources" au travers d'un "propagateur"
265!   (Green)
266! Elle montre aussi beaucoup plus clairement pourquoi on n'a pas ?? se souccier (trop)
267!   de la possible nulit?? du d??nominateur
268
269
270   do ilay=nlay,1,-1
271      do ig=1,ngrid
272         if (ilay.le.lmax(ig).and.lmax(ig)>1) then
273            edn(ig,ilay)=fact_thermals_down*dup(ig,ilay)
274            ddn(ig,ilay)=fact_thermals_down*eup(ig,ilay)
275            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
276            thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
277         endif
278      enddo
279   enddo
280
281   ! Suite du travail :
282   ! Ecrire la conservervation de theta_l dans le panache descendant
283   ! Eventuellement faire la transformation theta_l -> theta_v
284   ! Si l'air est sec (et qu'on oublie le c??t?? theta_v) on peut
285   ! se contenter de conserver theta.
286   !
287   ! Connaissant thetadn, on peut calculer la flotabilit??.
288   ! Connaissant la flotabilit??, on peut calculer w de proche en proche
289   ! On peut calculer le detrainement de facon ?? garder alpha*rho = cste
290   ! On en d??duit l'entrainement lat??ral
291   ! C'est le mod??le des mini-projets.
292
293!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
294! Initialisations :
295!------------------
296
297
298!
299 RETURN
300   END
Note: See TracBrowser for help on using the repository browser.