source: LMDZ6/branches/contrails/libf/phylmd/lmdz_thermcell_down.f90 @ 5450

Last change on this file since 5450 was 5390, checked in by yann meurdesoif, 3 weeks ago
  • Remove UTF8 character that inihibit fortran parsing with GPU morphosis
  • Add missing END SUBROUTINE instead of simple END, that inhibit correct parsing with regulat expression parser (quick and dirty parsing)

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