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

Last change on this file since 5135 was 5119, checked in by abarral, 16 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.