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

Last change on this file since 2176 was 2143, checked in by aboissinot, 6 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: 7.8 KB
Line 
1!
2!
3!
4SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2,             &
5                            zlev,zmin,zmix,zmax,wmax,f_star)
6     
7     
8!===============================================================================
9!  Purpose: calcul des caracteristiques du thermique: zmax,wmax,zmix
10!===============================================================================
11     
12      IMPLICIT NONE
13     
14     
15!===============================================================================
16! Declaration
17!===============================================================================
18     
19!     Inputs:
20!     -------
21     
22      INTEGER ngrid, nlay
23      INTEGER lmin(ngrid)
24     
25      REAL zlev(ngrid,nlay+1)
26      REAL f_star(ngrid,nlay+1)
27     
28!     Outputs:
29!     --------
30     
31      INTEGER lmix(ngrid)
32      INTEGER lmax(ngrid)
33     
34      REAL linter(ngrid)
35      REAL zmin(ngrid)                          ! Altitude of the plume bottom
36      REAL zmax(ngrid)                          ! Altitude of the plume top
37      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
38      REAL wmax(ngrid)                          ! Maximal vertical speed
39      REAL zw2(ngrid,nlay+1)                    ! Square vertical speed (becomes its square root)
40     
41!     Local:
42!     ------
43     
44      INTEGER ig, l
45     
46      REAL zlevinter(ngrid)
47     
48!===============================================================================
49! Initialization
50!===============================================================================
51     
52      DO ig=1,ngrid
53         lmax(ig) = lmin(ig)
54         lmix(ig) = lmin(ig)
55      ENDDO
56     
57      DO ig=1,ngrid
58         wmax(ig) = 0.
59      ENDDO
60     
61      DO ig=1,ngrid
62         zmax(ig) = 0.
63         zlevinter(ig) = zlev(ig,1)
64      ENDDO
65     
66!===============================================================================
67! Thermal plume height
68!===============================================================================
69     
70!-------------------------------------------------------------------------------
71! Calcul de zmin
72!-------------------------------------------------------------------------------
73     
74      DO l=1,nlay
75         DO ig=1,ngrid
76            zmin(ig) = zlev(ig,lmin(ig))
77         ENDDO
78      ENDDO
79     
80!-------------------------------------------------------------------------------
81! Calcul de lmax
82!-------------------------------------------------------------------------------
83     
84      DO ig=1,ngrid
85         DO l=nlay,lmin(ig)+1,-1
86            IF ((zw2(ig,l) <= 0.).or.(f_star(ig,l) <= 0.)) THEN
87               lmax(ig) = l - 1
88            ENDIF
89         ENDDO
90      ENDDO
91     
92!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
93! AB: Problematic case where thermal plume reaches the top of the model...
94!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
95      DO ig=1,ngrid
96         IF (zw2(ig,nlay) > 0.) THEN
97            print *, 'WARNING: a thermal plume reaches the highest layer!'
98            print *, 'ig,l', ig, nlay
99            print *, 'zw2', zw2(ig,nlay)
100            lmax(ig) = nlay
101         ENDIF
102      ENDDO
103     
104!-------------------------------------------------------------------------------
105! Calcul de zmax continu via le linter
106!-------------------------------------------------------------------------------
107     
108!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
109! AB : lmin=lmax means the plume is not active and then zw2=0 at each level.
110!      Otherwise we have lmax>lmin.
111!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
112      DO ig=1,ngrid
113         l = lmax(ig)
114         IF (l == lmin(ig)) THEN
115            linter(ig) = l
116         ELSE
117            linter(ig) = l - zw2(ig,l) / (zw2(ig,l+1) - zw2(ig,l))
118         ENDIF
119      ENDDO
120     
121      DO ig=1,ngrid
122         zlevinter(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig))          &
123         &             * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig)))
124         zmax(ig) = max(zmax(ig), zlevinter(ig))
125      ENDDO
126     
127!===============================================================================
128! Thermal plume maximal speed and inversion layer
129!===============================================================================
130     
131!-------------------------------------------------------------------------------
132! Calcul de lmix et wmax
133!-------------------------------------------------------------------------------
134     
135      DO l=1,nlay
136         DO ig=1,ngrid
137            IF((l <= lmax(ig)).and.(l > lmin(ig))) THEN
138               IF (zw2(ig,l).lt.0.) THEN
139                  print *, 'ERROR: zw2 has negative value(s)!'
140                  print *, 'ig,l', ig, l
141                  print *, 'zw2', zw2(ig,l)
142                  CALL abort
143               ENDIF
144               
145!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
146! AB : WARNING zw2 becomes its square root!
147!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
148               zw2(ig,l) = sqrt(zw2(ig,l))
149               
150               IF (zw2(ig,l) > wmax(ig)) THEN
151                  wmax(ig)  = zw2(ig,l)
152                  lmix(ig)  = l
153               ENDIF
154            ELSE
155               zw2(ig,l) = 0.
156            ENDIF
157         ENDDO
158      ENDDO
159     
160!-------------------------------------------------------------------------------
161! Calcul de zmix continu (profil parabolique des vitesses)
162!-------------------------------------------------------------------------------
163     
164      DO ig=1,ngrid
165         IF (lmix(ig) > lmin(ig)) THEN
166            IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))                        &
167            &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))             &
168            &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                   &
169            &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))) > 1e-10)   &
170            &        THEN
171               zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))              &
172               &        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)    &
173               &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                &
174               &        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))   &
175               &        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))           &
176               &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))          &
177               &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                &
178               &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
179            ELSE
180               zmix(ig) = zlev(ig,lmix(ig))
181               print *, 'WARNING: problematic zmix value!'
182            ENDIF
183         ELSE
184            zmix(ig) = 0.
185         ENDIF
186      ENDDO
187     
188!===============================================================================
189! Check consistence between zmax and zmix
190!===============================================================================
191     
192!-------------------------------------------------------------------------------
193!
194!-------------------------------------------------------------------------------
195     
196      DO ig=1,ngrid
197         IF ((zmax(ig)-zmix(ig)) < 0.) THEN
198            print *, 'WARNING: we have zmix > zmax!'
199            print *, 'zmax,zmix_old,zmix_new', zmax(ig), zmix(ig), 0.9 * zmax(ig)
200            zmix(ig) = 0.9 * zmax(ig)
201         ENDIF
202      ENDDO
203     
204!-------------------------------------------------------------------------------
205! Calcul du nouveau lmix correspondant
206!-------------------------------------------------------------------------------
207     
208      DO ig=1,ngrid
209         DO l=1,nlay
210            IF ((zmix(ig) >= zlev(ig,l)).and.(zmix(ig) < zlev(ig,l+1))) THEN
211               lmix(ig) = l
212            ENDIF
213         ENDDO
214      ENDDO
215     
216     
217RETURN
218END
Note: See TracBrowser for help on using the repository browser.