source: LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_dq.F90 @ 5119

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

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way
  • 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 lmdz_print_control, 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      CALL thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, &
48              masse, q, dq, qa, lev_out)
49
50    else
51
52
53      ! Calcul du critere CFL pour l'advection dans la subsidence
54      cfl = 0.
55      do k = 1, nlay
56        do ig = 1, ngrid
57          zzm = masse(ig, k) / ptimestep
58          cfl = max(cfl, fm(ig, k) / zzm)
59          IF (entr(ig, k)>zzm) THEN
60            PRINT*, 'entr*dt>m,1', k, entr(ig, k) * ptimestep, masse(ig, k)
61            abort_message = 'entr dt > m, 1st'
62            CALL abort_physic (modname, abort_message, 1)
63          endif
64        enddo
65      enddo
66
67      qold = q
68
69      IF (prt_level>=1) PRINT*, 'Q2 THERMCEL_DQ 0'
70
71      !   calcul du detrainement
72      do k = 1, nlay
73        do ig = 1, ngrid
74          detr(ig, k) = fm(ig, k) - fm(ig, k + 1) + entr(ig, k)
75          !           PRINT*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
76          !test
77          IF (detr(ig, k)<0.) THEN
78            entr(ig, k) = entr(ig, k) - detr(ig, k)
79            detr(ig, k) = 0.
80            !               PRINT*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
81            !     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
82          endif
83          IF (fm(ig, k + 1)<0.) THEN
84            !               PRINT*,'fm2<0!!!'
85          endif
86          IF (entr(ig, k)<0.) THEN
87            !               PRINT*,'entr2<0!!!'
88          endif
89        enddo
90      enddo
91
92      ! Computation of tracer concentrations in the ascending plume
93      do ig = 1, ngrid
94        qa(ig, 1) = q(ig, 1)
95      enddo
96
97      do k = 2, nlay
98        do ig = 1, ngrid
99          IF ((fm(ig, k + 1) + detr(ig, k)) * ptimestep>  &
100                  1.e-5 * masse(ig, k)) THEN
101            qa(ig, k) = (fm(ig, k) * qa(ig, k - 1) + entr(ig, k) * q(ig, k))  &
102                    / (fm(ig, k + 1) + detr(ig, k))
103          else
104            qa(ig, k) = q(ig, k)
105          endif
106          IF (qa(ig, k)<0.) THEN
107            !               PRINT*,'qa<0!!!'
108          endif
109          IF (q(ig, k)<0.) THEN
110            !               PRINT*,'q<0!!!'
111          endif
112        enddo
113      enddo
114
115      ! Plume vertical flux
116      do k = 2, nlay - 1
117        fqa(:, k) = fm(:, k) * qa(:, k - 1)
118      enddo
119      fqa(:, 1) = 0. ; fqa(:, nlay) = 0.
120
121
122      ! Trace species evolution
123      IF (impl==0) THEN
124        do k = 1, nlay - 1
125          q(:, k) = q(:, k) + (fqa(:, k) - fqa(:, k + 1) - fm(:, k) * q(:, k) + fm(:, k + 1) * q(:, k + 1)) &
126                  * ptimestep / masse(:, k)
127        enddo
128      else
129        do k = nlay - 1, 1, -1
130          ! FH debut de modif : le calcul ci dessous modifiait numériquement
131          ! la concentration quand le flux de masse etait nul car on divisait
132          ! puis multipliait par masse/ptimestep.
133          !        q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
134          !    &               /(fm(:,k)+masse(:,k)/ptimestep)
135          q(:, k) = (q(:, k) + ptimestep / masse(:, k) * (fqa(:, k) - fqa(:, k + 1) + fm(:, k + 1) * q(:, k + 1))) &
136                  / (1. + fm(:, k) * ptimestep / masse(:, k))
137          ! FH fin de modif.
138        enddo
139      endif
140
141      ! Tendencies
142      do k = 1, nlay
143        do ig = 1, ngrid
144          dq(ig, k) = (q(ig, k) - qold(ig, k)) / ptimestep
145          q(ig, k) = qold(ig, k)
146        enddo
147      enddo
148
149    END IF ! impl=-1
150    RETURN
151  end
152
153
154  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
155  ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
156  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
157
158  SUBROUTINE thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, &
159          masse, q, dq, qa, lev_out)
160    USE lmdz_print_control, ONLY: prt_level
161    USE lmdz_abort_physic, ONLY: abort_physic
162    IMPLICIT NONE
163
164    !=======================================================================
165
166    !   Calcul du transport verticale dans la couche limite en presence
167    !   de "thermiques" explicitement representes
168    !   calcul du dq/dt une fois qu'on connait les ascendances
169
170    !=======================================================================
171
172    INTEGER ngrid, nlay, impl
173
174    REAL ptimestep
175    REAL masse(ngrid, nlay), fm(ngrid, nlay + 1)
176    REAL entr(ngrid, nlay)
177    REAL q(ngrid, nlay)
178    REAL dq(ngrid, nlay)
179    INTEGER lev_out                           ! niveau pour les print
180
181    REAL qa(ngrid, nlay), detr(ngrid, nlay), wqd(ngrid, nlay + 1)
182
183    REAL zzm
184
185    INTEGER ig, k
186    REAL cfl
187
188    REAL qold(ngrid, nlay)
189    REAL ztimestep
190    INTEGER niter, iter
191    CHARACTER (LEN = 20) :: modname = 'thermcell_dq'
192    CHARACTER (LEN = 80) :: abort_message
193
194
195
196    ! Calcul du critere CFL pour l'advection dans la subsidence
197    cfl = 0.
198    do k = 1, nlay
199      do ig = 1, ngrid
200        zzm = masse(ig, k) / ptimestep
201        cfl = max(cfl, fm(ig, k) / zzm)
202        IF (entr(ig, k)>zzm) THEN
203          PRINT*, 'entr*dt>m,2', k, entr(ig, k) * ptimestep, masse(ig, k)
204          abort_message = 'entr dt > m, 2nd'
205          CALL abort_physic (modname, abort_message, 1)
206        endif
207      enddo
208    enddo
209
210    !IM 090508     PRINT*,'CFL CFL CFL CFL ',cfl
211
212#undef CFL
213#ifdef CFL
214! On subdivise le calcul en niter pas de temps.
215      niter=int(cfl)+1
216#else
217    niter = 1
218#endif
219
220    ztimestep = ptimestep / niter
221    qold = q
222
223    DO iter = 1, niter
224      IF (prt_level>=1) PRINT*, 'Q2 THERMCEL_DQ 0'
225
226      !   calcul du detrainement
227      do k = 1, nlay
228        do ig = 1, ngrid
229          detr(ig, k) = fm(ig, k) - fm(ig, k + 1) + entr(ig, k)
230          !           PRINT*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
231          !test
232          IF (detr(ig, k)<0.) THEN
233            entr(ig, k) = entr(ig, k) - detr(ig, k)
234            detr(ig, k) = 0.
235            !               PRINT*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
236            !     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
237          endif
238          IF (fm(ig, k + 1)<0.) THEN
239            !               PRINT*,'fm2<0!!!'
240          endif
241          IF (entr(ig, k)<0.) THEN
242            !               PRINT*,'entr2<0!!!'
243          endif
244        enddo
245      enddo
246
247      !   calcul de la valeur dans les ascendances
248      do ig = 1, ngrid
249        qa(ig, 1) = q(ig, 1)
250      enddo
251
252      do k = 2, nlay
253        do ig = 1, ngrid
254          IF ((fm(ig, k + 1) + detr(ig, k)) * ztimestep>  &
255                  1.e-5 * masse(ig, k)) THEN
256            qa(ig, k) = (fm(ig, k) * qa(ig, k - 1) + entr(ig, k) * q(ig, k))  &
257                    / (fm(ig, k + 1) + detr(ig, k))
258          else
259            qa(ig, k) = q(ig, k)
260          endif
261          IF (qa(ig, k)<0.) THEN
262            !               PRINT*,'qa<0!!!'
263          endif
264          IF (q(ig, k)<0.) THEN
265            !               PRINT*,'q<0!!!'
266          endif
267        enddo
268      enddo
269
270      ! Calcul du flux subsident
271
272      do k = 2, nlay
273        do ig = 1, ngrid
274#undef centre
275#ifdef centre
276             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
277#else
278
279#define CFL_plus_grand_que_un
280#ifdef CFL_plus_grand_que_un
281          ! Schema avec advection sur plus qu'une maille.
282          zzm = masse(ig, k) / ztimestep
283          IF (fm(ig, k)>zzm) THEN
284            wqd(ig, k) = zzm * q(ig, k) + (fm(ig, k) - zzm) * q(ig, k + 1)
285          else
286            wqd(ig, k) = fm(ig, k) * q(ig, k)
287          endif
288#else
289            wqd(ig,k)=fm(ig,k)*q(ig,k)
290#endif
291#endif
292
293          IF (wqd(ig, k)<0.) THEN
294            !               PRINT*,'wqd<0!!!'
295          endif
296        enddo
297      enddo
298      do ig = 1, ngrid
299        wqd(ig, 1) = 0.
300        wqd(ig, nlay + 1) = 0.
301      enddo
302
303
304      ! Calcul des tendances
305      do k = 1, nlay
306        do ig = 1, ngrid
307          q(ig, k) = q(ig, k) + (detr(ig, k) * qa(ig, k) - entr(ig, k) * q(ig, k)  &
308                  - wqd(ig, k) + wqd(ig, k + 1))  &
309                  * ztimestep / masse(ig, k)
310          !            if (dq(ig,k).lt.0.) THEN
311          !               PRINT*,'dq<0!!!'
312          !            endif
313        enddo
314      enddo
315
316    END DO
317
318
319    ! Calcul des tendances
320    do k = 1, nlay
321      do ig = 1, ngrid
322        dq(ig, k) = (q(ig, k) - qold(ig, k)) / ptimestep
323        q(ig, k) = qold(ig, k)
324      enddo
325    enddo
326
327    RETURN
328  END
329END MODULE lmdz_thermcell_dq
Note: See TracBrowser for help on using the repository browser.