source: trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90 @ 2176

Last change on this file since 2176 was 2144, checked in by aboissinot, 5 years ago

Thermal plume model parameters are now set in callphys.def and (re)add flag to choose between thermcell_dv2 or thermcell_dq to mix horizontal momentum.
File size: 5.2 KB
Line 
1!
2!
3!
4SUBROUTINE thermcell_dq(ngrid,nlay,ptimestep,fm,entr,detr,masse,              &
5                        q,dq,qa,lmin)
6     
7     
8!===============================================================================
9!  Purpose: Calcul du transport verticale dans la couche limite en presence de
10!           "thermiques" explicitement representes
11!           Calcul du dq/dt une fois qu'on connait les ascendances
12! 
13!  Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
14!  Introduction of an implicit computation of vertical advection in the environ-
15!     ment of thermal plumes in thermcell_dq
16! 
17!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
18!     dqimpl = 1 : implicit scheme
19!     dqimpl = 0 : explicit scheme
20! 
21!===============================================================================
22     
23      USE print_control_mod, ONLY: prt_level
24      USe thermcell_mod, ONLY: dqimpl
25     
26      IMPLICIT NONE
27     
28     
29!===============================================================================
30! Declaration
31!===============================================================================
32     
33!     Inputs:
34!     -------
35     
36      INTEGER ngrid, nlay
37      INTEGER lmin(ngrid)
38     
39      REAL ptimestep
40      REAL masse(ngrid,nlay)
41      REAL fm(ngrid,nlay+1)
42      REAL entr(ngrid,nlay)
43      REAL detr(ngrid,nlay)
44      REAL q(ngrid,nlay)
45     
46!     Outputs:
47!     --------
48     
49      REAL dq(ngrid,nlay)
50      REAL qa(ngrid,nlay)
51     
52!     Local:
53!     ------
54     
55      INTEGER ig, l, k
56      INTEGER niter, iter
57     
58      REAL cfl
59      REAL qold(ngrid,nlay)
60      REAL fqa(ngrid,nlay+1)
61      REAL zzm
62     
63!===============================================================================
64! Initialization
65!===============================================================================
66     
67      qold(:,:) = q(:,:)
68     
69!===============================================================================
70! Tracer variation computation
71!===============================================================================
72     
73!-------------------------------------------------------------------------------
74! CFL criterion computation for advection in downdraft
75!-------------------------------------------------------------------------------
76     
77      cfl = 0.
78     
79      DO l=1,nlay
80         DO ig=1,ngrid
81            zzm = masse(ig,l) / ptimestep
82            cfl = max(cfl, fm(ig,l) / zzm)
83           
84            IF (entr(ig,l).gt.zzm) THEN
85               print *, 'ERROR: entrainment is greater than the layer mass!'
86               print *, 'ig,l,entr', ig, l, entr(ig,l)
87               print *, 'lmin', lmin(ig)
88               print *, '-------------------------------'
89               print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l)
90               print *, '-------------------------------'
91               DO k=nlay,1,-1
92                  print *, 'fm ', fm(ig,k+1)
93                  print *, 'entr,detr', entr(ig,k), detr(ig,k)
94               ENDDO
95               print *, 'fm ', fm(ig,1)
96               print *, '-------------------------------'
97               CALL abort
98            ENDIF
99         ENDDO
100      ENDDO
101     
102!-------------------------------------------------------------------------------
103! Computation of tracer concentrations in the ascending plume
104!-------------------------------------------------------------------------------
105     
106      DO ig=1,ngrid
107         DO l=1,lmin(ig)
108            qa(ig,l) = q(ig,l)
109         ENDDO
110      ENDDO
111     
112      DO ig=1,ngrid
113         DO l=lmin(ig)+1,nlay
114            IF ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-6*masse(ig,l)) THEN
115               qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l))      &
116               &        / (fm(ig,l+1) + detr(ig,l))
117            ELSE
118               qa(ig,l) = q(ig,l)
119            ENDIF
120         ENDDO
121      ENDDO
122     
123!-------------------------------------------------------------------------------
124! Plume vertical flux of tracer
125!-------------------------------------------------------------------------------
126     
127      DO l=2,nlay-1
128         fqa(:,l) = fm(:,l) * qa(:,l-1)
129      ENDDO
130     
131      fqa(:,1) = 0.
132      fqa(:,nlay) = 0.
133     
134!-------------------------------------------------------------------------------
135! Trace species evolution
136!-------------------------------------------------------------------------------
137     
138      IF (dqimpl) THEN
139         DO l=nlay-1,1,-1
140            q(:,l) = ( q(:,l) + ptimestep / masse(:,l)                        &
141            &      * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) )       &
142            &      / ( 1. + fm(:,l) * ptimestep / masse(:,l) )
143         ENDDO
144      ELSE
145         DO l=1,nlay-1
146            q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l)       &
147            &      + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l)
148         ENDDO
149      ENDIF
150     
151!===============================================================================
152! Tendencies
153!===============================================================================
154     
155      DO l=1,nlay
156         DO ig=1,ngrid
157            dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
158            q(ig,l) = qold(ig,l)
159         ENDDO
160      ENDDO
161     
162     
163RETURN
164END
Note: See TracBrowser for help on using the repository browser.