source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_dq.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
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.6 KB
Line 
1MODULE lmdz_thermcell_dq
2CONTAINS
3
4  SUBROUTINE thermcell_dq(ngrid, nlay, impl, ptimestep, fm, entr, &
5          masse, q, dq, qa, lev_out)
6    USE print_control_mod, ONLY: prt_level
7    USE lmdz_abort_physic, ONLY: abort_physic
8
9    implicit none
10
11    !=======================================================================
12
13    !   Calcul du transport verticale dans la couche limite en presence
14    !   de "thermiques" explicitement representes
15    !   calcul du dq/dt une fois qu'on connait les ascendances
16
17    ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
18    !  Introduction of an implicit computation of vertical advection in
19    !  the environment of thermal plumes in thermcell_dq
20    !  impl =     0 : explicit, 1 : implicit, -1 : old version
21
22    !=======================================================================
23
24    ! arguments
25    integer, intent(in) :: ngrid, nlay, impl
26    real, intent(in) :: ptimestep
27    real, intent(in), dimension(ngrid, nlay) :: masse
28    real, intent(inout), dimension(ngrid, nlay) :: entr, q
29    real, intent(in), dimension(ngrid, nlay + 1) :: fm
30    real, intent(out), dimension(ngrid, nlay) :: dq, qa
31    integer, intent(in) :: lev_out                           ! niveau pour les print
32
33    ! Local
34    real, dimension(ngrid, nlay) :: detr, qold
35    real, dimension(ngrid, nlay + 1) :: wqd, fqa
36    real zzm
37    integer ig, k
38    real cfl
39
40    integer niter, iter
41    CHARACTER (LEN = 20) :: modname = 'thermcell_dq'
42    CHARACTER (LEN = 80) :: abort_message
43
44
45    ! Old explicite scheme
46    IF (impl<=-1) then
47
48      CALL thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, &
49              masse, q, dq, qa, lev_out)
50
51    else
52
53
54      ! Calcul du critere CFL pour l'advection dans la subsidence
55      cfl = 0.
56      do k = 1, nlay
57        do ig = 1, ngrid
58          zzm = masse(ig, k) / ptimestep
59          cfl = max(cfl, fm(ig, k) / zzm)
60          if (entr(ig, k)>zzm) then
61            PRINT*, 'entr*dt>m,1', k, entr(ig, k) * ptimestep, masse(ig, k)
62            abort_message = 'entr dt > m, 1st'
63            CALL abort_physic (modname, abort_message, 1)
64          endif
65        enddo
66      enddo
67
68      qold = q
69
70      if (prt_level>=1) PRINT*, 'Q2 THERMCEL_DQ 0'
71
72      !   calcul du detrainement
73      do k = 1, nlay
74        do ig = 1, ngrid
75          detr(ig, k) = fm(ig, k) - fm(ig, k + 1) + entr(ig, k)
76          !           PRINT*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
77          !test
78          if (detr(ig, k)<0.) then
79            entr(ig, k) = entr(ig, k) - detr(ig, k)
80            detr(ig, k) = 0.
81            !               PRINT*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
82            !     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
83          endif
84          if (fm(ig, k + 1)<0.) then
85            !               PRINT*,'fm2<0!!!'
86          endif
87          if (entr(ig, k)<0.) then
88            !               PRINT*,'entr2<0!!!'
89          endif
90        enddo
91      enddo
92
93      ! Computation of tracer concentrations in the ascending plume
94      do ig = 1, ngrid
95        qa(ig, 1) = q(ig, 1)
96      enddo
97
98      do k = 2, nlay
99        do ig = 1, ngrid
100          if ((fm(ig, k + 1) + detr(ig, k)) * ptimestep>  &
101                  1.e-5 * masse(ig, k)) then
102            qa(ig, k) = (fm(ig, k) * qa(ig, k - 1) + entr(ig, k) * q(ig, k))  &
103                    / (fm(ig, k + 1) + detr(ig, k))
104          else
105            qa(ig, k) = q(ig, k)
106          endif
107          if (qa(ig, k)<0.) then
108            !               PRINT*,'qa<0!!!'
109          endif
110          if (q(ig, k)<0.) then
111            !               PRINT*,'q<0!!!'
112          endif
113        enddo
114      enddo
115
116      ! Plume vertical flux
117      do k = 2, nlay - 1
118        fqa(:, k) = fm(:, k) * qa(:, k - 1)
119      enddo
120      fqa(:, 1) = 0. ; fqa(:, nlay) = 0.
121
122
123      ! Trace species evolution
124      if (impl==0) then
125        do k = 1, nlay - 1
126          q(:, k) = q(:, k) + (fqa(:, k) - fqa(:, k + 1) - fm(:, k) * q(:, k) + fm(:, k + 1) * q(:, k + 1)) &
127                  * ptimestep / masse(:, k)
128        enddo
129      else
130        do k = nlay - 1, 1, -1
131          ! FH debut de modif : le calcul ci dessous modifiait numériquement
132          ! la concentration quand le flux de masse etait nul car on divisait
133          ! puis multipliait par masse/ptimestep.
134          !        q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
135          !    &               /(fm(:,k)+masse(:,k)/ptimestep)
136          q(:, k) = (q(:, k) + ptimestep / masse(:, k) * (fqa(:, k) - fqa(:, k + 1) + fm(:, k + 1) * q(:, k + 1))) &
137                  / (1. + fm(:, k) * ptimestep / masse(:, k))
138          ! FH fin de modif.
139        enddo
140      endif
141
142      ! Tendencies
143      do k = 1, nlay
144        do ig = 1, ngrid
145          dq(ig, k) = (q(ig, k) - qold(ig, k)) / ptimestep
146          q(ig, k) = qold(ig, k)
147        enddo
148      enddo
149
150    END IF ! impl=-1
151    RETURN
152  end
153
154
155  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
156  ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
157  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158
159  SUBROUTINE thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, &
160          masse, q, dq, qa, lev_out)
161    USE print_control_mod, ONLY: prt_level
162    USE lmdz_abort_physic, ONLY: abort_physic
163    implicit none
164
165    !=======================================================================
166
167    !   Calcul du transport verticale dans la couche limite en presence
168    !   de "thermiques" explicitement representes
169    !   calcul du dq/dt une fois qu'on connait les ascendances
170
171    !=======================================================================
172
173    integer ngrid, nlay, impl
174
175    real ptimestep
176    real masse(ngrid, nlay), fm(ngrid, nlay + 1)
177    real entr(ngrid, nlay)
178    real q(ngrid, nlay)
179    real dq(ngrid, nlay)
180    integer lev_out                           ! niveau pour les print
181
182    real qa(ngrid, nlay), detr(ngrid, nlay), wqd(ngrid, nlay + 1)
183
184    real zzm
185
186    integer ig, k
187    real cfl
188
189    real qold(ngrid, nlay)
190    real ztimestep
191    integer niter, iter
192    CHARACTER (LEN = 20) :: modname = 'thermcell_dq'
193    CHARACTER (LEN = 80) :: abort_message
194
195
196
197    ! Calcul du critere CFL pour l'advection dans la subsidence
198    cfl = 0.
199    do k = 1, nlay
200      do ig = 1, ngrid
201        zzm = masse(ig, k) / ptimestep
202        cfl = max(cfl, fm(ig, k) / zzm)
203        if (entr(ig, k)>zzm) then
204          PRINT*, 'entr*dt>m,2', k, entr(ig, k) * ptimestep, masse(ig, k)
205          abort_message = 'entr dt > m, 2nd'
206          CALL abort_physic (modname, abort_message, 1)
207        endif
208      enddo
209    enddo
210
211    !IM 090508     PRINT*,'CFL CFL CFL CFL ',cfl
212
213#undef CFL
214#ifdef CFL
215! On subdivise le calcul en niter pas de temps.
216      niter=int(cfl)+1
217#else
218    niter = 1
219#endif
220
221    ztimestep = ptimestep / niter
222    qold = q
223
224    DO iter = 1, niter
225      if (prt_level>=1) PRINT*, 'Q2 THERMCEL_DQ 0'
226
227      !   calcul du detrainement
228      do k = 1, nlay
229        do ig = 1, ngrid
230          detr(ig, k) = fm(ig, k) - fm(ig, k + 1) + entr(ig, k)
231          !           PRINT*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
232          !test
233          if (detr(ig, k)<0.) then
234            entr(ig, k) = entr(ig, k) - detr(ig, k)
235            detr(ig, k) = 0.
236            !               PRINT*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
237            !     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
238          endif
239          if (fm(ig, k + 1)<0.) then
240            !               PRINT*,'fm2<0!!!'
241          endif
242          if (entr(ig, k)<0.) then
243            !               PRINT*,'entr2<0!!!'
244          endif
245        enddo
246      enddo
247
248      !   calcul de la valeur dans les ascendances
249      do ig = 1, ngrid
250        qa(ig, 1) = q(ig, 1)
251      enddo
252
253      do k = 2, nlay
254        do ig = 1, ngrid
255          if ((fm(ig, k + 1) + detr(ig, k)) * ztimestep>  &
256                  1.e-5 * masse(ig, k)) then
257            qa(ig, k) = (fm(ig, k) * qa(ig, k - 1) + entr(ig, k) * q(ig, k))  &
258                    / (fm(ig, k + 1) + detr(ig, k))
259          else
260            qa(ig, k) = q(ig, k)
261          endif
262          if (qa(ig, k)<0.) then
263            !               PRINT*,'qa<0!!!'
264          endif
265          if (q(ig, k)<0.) then
266            !               PRINT*,'q<0!!!'
267          endif
268        enddo
269      enddo
270
271      ! Calcul du flux subsident
272
273      do k = 2, nlay
274        do ig = 1, ngrid
275#undef centre
276#ifdef centre
277             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
278#else
279
280#define CFL_plus_grand_que_un
281#ifdef CFL_plus_grand_que_un
282          ! Schema avec advection sur plus qu'une maille.
283          zzm = masse(ig, k) / ztimestep
284          if (fm(ig, k)>zzm) then
285            wqd(ig, k) = zzm * q(ig, k) + (fm(ig, k) - zzm) * q(ig, k + 1)
286          else
287            wqd(ig, k) = fm(ig, k) * q(ig, k)
288          endif
289#else
290            wqd(ig,k)=fm(ig,k)*q(ig,k)
291#endif
292#endif
293
294          if (wqd(ig, k)<0.) then
295            !               PRINT*,'wqd<0!!!'
296          endif
297        enddo
298      enddo
299      do ig = 1, ngrid
300        wqd(ig, 1) = 0.
301        wqd(ig, nlay + 1) = 0.
302      enddo
303
304
305      ! Calcul des tendances
306      do k = 1, nlay
307        do ig = 1, ngrid
308          q(ig, k) = q(ig, k) + (detr(ig, k) * qa(ig, k) - entr(ig, k) * q(ig, k)  &
309                  - wqd(ig, k) + wqd(ig, k + 1))  &
310                  * ztimestep / masse(ig, k)
311          !            if (dq(ig,k).lt.0.) then
312          !               PRINT*,'dq<0!!!'
313          !            endif
314        enddo
315      enddo
316
317    END DO
318
319
320    ! Calcul des tendances
321    do k = 1, nlay
322      do ig = 1, ngrid
323        dq(ig, k) = (q(ig, k) - qold(ig, k)) / ptimestep
324        q(ig, k) = qold(ig, k)
325      enddo
326    enddo
327
328    return
329  end
330END MODULE lmdz_thermcell_dq
Note: See TracBrowser for help on using the repository browser.