source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_down.F90 @ 5111

Last change on this file since 5111 was 5111, checked in by abarral, 4 months ago

Put abort_physic into a module
Remove -g option from makelmdz_fcm, since that option is linked to a header file that isn't included anywhere.
(lint) light lint on traversed files
File size: 11.5 KB
Line 
1MODULE lmdz_thermcell_down
2CONTAINS
3
4  SUBROUTINE thermcell_updown_dq(ngrid, nlay, ptimestep, lmax, eup, dup, edn, ddn, masse, trac, dtrac)
5
6    USE 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    USE lmdz_abort_physic, ONLY: abort_physic
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<=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<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) < 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) < 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)) < 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
216  END
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    USE lmdz_thermcell_ini, ONLY: prt_level, RLvCp, RKAPPA, RETV, fact_thermals_down
228    IMPLICIT NONE
229
230    ! arguments
231
232    integer, intent(in) :: ngrid, nlay
233    real, intent(in), dimension(ngrid, nlay) :: po, pt, pu, pv, pplay, eup, dup
234    real, intent(in), dimension(ngrid, nlay) :: theta
235    real, intent(in), dimension(ngrid, nlay + 1) :: pplev, fup
236    integer, intent(in), dimension(ngrid) :: lmax
237
238
239
240    ! Local
241
242    real, dimension(ngrid, nlay) :: edn, ddn, thetad
243    real, dimension(ngrid, nlay + 1) :: fdn
244
245    integer ig, ilay
246    real dqsat_dT
247    logical mask(ngrid, nlay)
248
249    edn(:, :) = 0.
250    ddn(:, :) = 0.
251    fdn(:, :) = 0.
252    thetad(:, :) = 0.
253
254    ! lmax : indice tel que fu(kmax+1)=0
255
256    ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) )
257
258    ! FH MODIFS APRES REUNIONS POUR COMMISSIONS
259    ! quelques erreurs de declaration
260    ! probleme si lmax=1 ce qui a l'air d'??tre le cas en d??but de simu. Devrait ??tre 0 ?
261    ! Remarques :
262    ! on pourrait ??crire la formule de thetad
263    !    www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
264    !    thetad(ig,ilay)= www * thetad(ig,ilay+1) + (1.-www) * theta(ig,ilay)
265    ! Elle a l'avantage de bien montr?? la conservation, l'id??e fondamentale dans le
266    !   transport qu'on ne fait que sommer des "sources" au travers d'un "propagateur"
267    !   (Green)
268    ! Elle montre aussi beaucoup plus clairement pourquoi on n'a pas ?? se souccier (trop)
269    !   de la possible nulit?? du d??nominateur
270
271    do ilay = nlay, 1, -1
272      do ig = 1, ngrid
273        if (ilay<=lmax(ig).and.lmax(ig)>1) then
274          edn(ig, ilay) = fact_thermals_down * dup(ig, ilay)
275          ddn(ig, ilay) = fact_thermals_down * eup(ig, ilay)
276          fdn(ig, ilay) = fdn(ig, ilay + 1) + edn(ig, ilay) - ddn(ig, ilay)
277          thetad(ig, ilay) = (fdn(ig, ilay + 1) * thetad(ig, ilay + 1) + edn(ig, ilay) * theta(ig, ilay)) / (fdn(ig, ilay) + ddn(ig, ilay))
278        endif
279      enddo
280    enddo
281
282    ! Suite du travail :
283    ! Ecrire la conservervation de theta_l dans le panache descendant
284    ! Eventuellement faire la transformation theta_l -> theta_v
285    ! Si l'air est sec (et qu'on oublie le c??t?? theta_v) on peut
286    ! se contenter de conserver theta.
287
288    ! Connaissant thetadn, on peut calculer la flotabilit??.
289    ! Connaissant la flotabilit??, on peut calculer w de proche en proche
290    ! On peut calculer le detrainement de facon ?? garder alpha*rho = cste
291    ! On en d??duit l'entrainement lat??ral
292    ! C'est le mod??le des mini-projets.
293
294    !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
295    ! Initialisations :
296    !------------------
297
298    RETURN
299  END
300END MODULE lmdz_thermcell_down
Note: See TracBrowser for help on using the repository browser.