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

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

Update the thermal plume model:

  • check formulae consistency between thermcell_plume and thercell_closure
  • compute correctly thermal plume height
  • fix alimentation computation in the first unstable layer
File size: 6.1 KB
Line 
1!
2!
3!
4SUBROUTINE thermcell_dv2(ngrid,nlay,ptimestep,fm,entr,detr,masse,fraca,       &
5                         zmax,zmin,u,v,du,dv,ua,va)
6     
7     
8!===============================================================================
9!
10!   Calcul du transport verticale dans la couche limite en presence
11!   de "thermiques" explicitement representes
12!   calcul du dq/dt une fois qu'on connait les ascendances
13!
14! Vectorisation, FH : 2010/03/08
15!
16!===============================================================================
17     
18      USE print_control_mod, ONLY: prt_level,lunout
19     
20      IMPLICIT none
21     
22     
23!===============================================================================
24! Declaration
25!===============================================================================
26     
27!     Inputs:
28!     -------
29     
30      INTEGER ngrid, nlay
31     
32      REAL ptimestep
33      REAL masse(ngrid,nlay)
34      REAL fm(ngrid,nlay+1)
35      REAL entr(ngrid,nlay)
36      REAL detr(ngrid,nlay)
37      REAL fraca(ngrid,nlay+1)
38      REAL zmax(ngrid)
39      REAL zmin(ngrid)
40      REAL u(ngrid,nlay)
41      REAL v(ngrid,nlay)
42     
43!     Outputs:
44!     --------
45     
46      REAL dua(ngrid,nlay)
47      REAL dva(ngrid,nlay)
48      REAL ua(ngrid,nlay)
49      REAL va(ngrid,nlay)
50      REAL du(ngrid,nlay)
51      REAL dv(ngrid,nlay)
52     
53!     Local:
54!     ------
55     
56      INTEGER ig, l
57      INTEGER iter
58     
59      REAL qa(ngrid,nlay)
60      REAL zf
61      REAL zf2
62      REAL wvd(ngrid,nlay+1)
63      REAL wud(ngrid,nlay+1)
64      REAL gamma0(ngrid,nlay+1)
65      REAL gamma(ngrid,nlay+1)
66      REAL ue(ngrid,nlay)
67      REAL ve(ngrid,nlay)
68      REAL plume_height(ngrid)
69     
70      LOGICAL ltherm(ngrid,nlay)
71     
72!===============================================================================
73! Initialization
74!===============================================================================
75     
76      DO ig=1,ngrid
77         plume_height(ig) = zmax(ig) - zmin(ig)
78      ENDDO
79     
80      DO ig=1,ngrid
81         ua(ig,1)=u(ig,1)
82         va(ig,1)=v(ig,1)
83         ue(ig,1)=u(ig,1)
84         ve(ig,1)=v(ig,1)
85      ENDDO
86     
87      gamma0(1:ngrid,1)=0.
88     
89      DO l=2,nlay
90         DO ig=1,ngrid
91            ltherm(ig,l) = (fm(ig,l+1) + detr(ig,l)) * ptimestep > 1.e-5 * masse(ig,l)
92            IF (ltherm(ig,l).and.(zmax(ig) > 0.)) THEN
93               gamma0(ig,l) = masse(ig,l) * 0.5 / zmax(ig)                    &
94               &            * sqrt(0.5 * (fraca(ig,l+1) + fraca(ig,l)))
95            ELSE
96               gamma0(ig,l) = 0.
97            ENDIF
98         ENDDO
99      ENDDO
100     
101      gamma(:,:) = 0.
102     
103!===============================================================================
104!
105!===============================================================================
106     
107      DO l=2,nlay
108         
109         DO ig=1,ngrid
110            IF (ltherm(ig,l)) THEN
111               dua(ig,l)=ua(ig,l-1)-u(ig,l-1)
112               dva(ig,l)=va(ig,l-1)-v(ig,l-1)
113            ELSE
114               ua(ig,l)=u(ig,l)
115               va(ig,l)=v(ig,l)
116               ue(ig,l)=u(ig,l)
117               ve(ig,l)=v(ig,l)
118            ENDIF
119         ENDDO
120         
121         DO iter=1,5
122            DO ig=1,ngrid
123! Calcul prenant en compte la fraction reelle
124               zf = 0.5 * (fraca(ig,l) + fraca(ig,l+1))
125               zf2 = 1./(1.-zf)
126! Calcul avec fraction infiniement petite
127!               zf = 0.
128!               zf2 = 1.
129               
130!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
131! FH: La première fois on multiplie le coefficient de freinage par le module du
132!     vent dans la couche en dessous. Mais pourquoi donc ?
133!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
134               IF (ltherm(ig,l)) THEN
135!   On choisit une relaxation lineaire.
136!                  gamma(ig,l) = gamma0(ig,l)
137!   On choisit une relaxation quadratique.
138                  gamma(ig,l) = gamma0(ig,l) * sqrt(dua(ig,l)**2 + dva(ig,l)**2)
139                 
140                  ua(ig,l) = (fm(ig,l) * ua(ig,l-1)                           &
141                  &        + (zf2 * entr(ig,l) + gamma(ig,l)) * u(ig,l))      &
142                  &        / (fm(ig,l+1) + detr(ig,l) + entr(ig,l) * zf * zf2 &
143                  &        + gamma(ig,l))
144                  va(ig,l) = (fm(ig,l) * va(ig,l-1)                           &
145                  &        + (zf2 * entr(ig,l) + gamma(ig,l)) * v(ig,l))      &
146                  &        / (fm(ig,l+1) + detr(ig,l) + entr(ig,l) * zf * zf2 &
147                  &        + gamma(ig,l))
148                  dua(ig,l) = ua(ig,l) - u(ig,l)
149                  dva(ig,l) = va(ig,l) - v(ig,l)
150                  ue(ig,l) = (u(ig,l) - zf * ua(ig,l)) * zf2
151                  ve(ig,l) = (v(ig,l) - zf * va(ig,l)) * zf2
152               ENDIF
153            ENDDO
154         ENDDO
155      ENDDO
156     
157!-------------------------------------------------------------------------------
158! Calcul du flux vertical de moment dans l'environnement
159!-------------------------------------------------------------------------------
160     
161      DO l=2,nlay
162         DO ig=1,ngrid
163            wud(ig,l) = fm(ig,l) * ue(ig,l)
164            wvd(ig,l) = fm(ig,l) * ve(ig,l)
165         ENDDO
166      ENDDO
167     
168      DO ig=1,ngrid
169         wud(ig,1) = 0.
170         wud(ig,nlay+1) = 0.
171         wvd(ig,1) = 0.
172         wvd(ig,nlay+1) = 0.
173      ENDDO
174     
175!===============================================================================
176! Tendencies
177!===============================================================================
178     
179      DO l=1,nlay
180         DO ig=1,ngrid
181            du(ig,l) = ((detr(ig,l) + gamma(ig,l)) * ua(ig,l)                 &
182            &        - (entr(ig,l) + gamma(ig,l)) * ue(ig,l)                  &
183            &        - wud(ig,l) + wud(ig,l+1))                               &
184            &        / masse(ig,l)
185            dv(ig,l) = ((detr(ig,l) + gamma(ig,l)) * va(ig,l)                 &
186            &        - (entr(ig,l) + gamma(ig,l)) * ve(ig,l)                  &
187            &        - wvd(ig,l) + wvd(ig,l+1))                               &
188            &        / masse(ig,l)
189         ENDDO
190      ENDDO
191     
192     
193RETURN
194END
Note: See TracBrowser for help on using the repository browser.