source: LMDZ6/trunk/libf/phylmd/lmdz_thermcell_down.f90 @ 5523

Last change on this file since 5523 was 5512, checked in by yann meurdesoif, 8 days ago

Implement GPU automatic port for :

  • Thermics
  • acama_gwd_rando
  • flott_gwd_rando

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