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

Last change on this file since 3436 was 2480, checked in by aboissinot, 4 years ago

Commit the last changes in the thermal plume model (more english comments).
File size: 7.4 KB
Line 
1!
2!
3!
4SUBROUTINE thermcell_height(ngrid,nlay,lmin,lmix,lmax,                        &
5                            zlev,zmin,zmix,zmax,zw2,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, INTENT(in) :: ngrid
23      INTEGER, INTENT(in) :: nlay
24      INTEGER, INTENT(in) :: lmin(ngrid)
25     
26      REAL, INTENT(in) :: zlev(ngrid,nlay+1)
27      REAL, INTENT(in) :: f_star(ngrid,nlay+1)
28     
29!     Outputs:
30!     --------
31     
32      INTEGER, INTENT(out) :: lmix(ngrid)
33      INTEGER, INTENT(out) :: lmax(ngrid)
34     
35      REAL, INTENT(out) :: zmin(ngrid)                ! Altitude of the plume bottom
36      REAL, INTENT(out) :: zmix(ngrid)                ! Altitude of maximal vertical speed
37      REAL, INTENT(out) :: zmax(ngrid)                ! Altitude of the plume top
38      REAL, INTENT(out) :: wmax(ngrid)                ! Maximal vertical speed
39     
40      REAL, INTENT(inout) :: zw2(ngrid,nlay+1)        ! Square vertical speed (becomes its square root)
41     
42!     Local:
43!     ------
44     
45      INTEGER ig, k
46     
47      REAL linter(ngrid)                              ! Level (continuous) of maximal vertical speed
48     
49!===============================================================================
50! Initialization
51!===============================================================================
52     
53      linter(:) = lmin(:)
54      lmix(:) = lmin(:)
55      lmax(:) = lmin(:)
56     
57      wmax(:) = 0.
58     
59      zmin(:) = 0.
60      zmix(:) = 0.
61      zmax(:) = 0.
62     
63!===============================================================================
64! Thermal plume height
65!===============================================================================
66     
67!-------------------------------------------------------------------------------
68! Calcul de zmin
69!-------------------------------------------------------------------------------
70     
71      DO ig=1,ngrid
72         zmin(ig) = zlev(ig,lmin(ig))
73      ENDDO
74     
75!-------------------------------------------------------------------------------
76! Calcul de lmax
77!-------------------------------------------------------------------------------
78     
79      DO ig=1,ngrid
80         DO k=nlay,lmin(ig)+1,-1
81            IF ((zw2(ig,k) <= 0.).or.(f_star(ig,k) <= 0.)) THEN
82               lmax(ig) = k - 1
83            ENDIF
84         ENDDO
85      ENDDO
86     
87!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
88! AB: Problematic case where thermal plume reaches the top of the model...
89!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
90      DO ig=1,ngrid
91         IF (zw2(ig,nlay) > 0.) THEN
92            print *, 'WARNING: a thermal plume reaches the highest layer!'
93            print *, 'ig,k', ig, nlay
94            print *, 'zw2', zw2(ig,nlay)
95            lmax(ig) = nlay
96         ENDIF
97      ENDDO
98     
99!-------------------------------------------------------------------------------
100! Calcul de zmax continu via linter
101!-------------------------------------------------------------------------------
102     
103!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
104! AB: lmin=lmax means the plume is not active and then zw2=0 at each level.
105!     Otherwise we have lmax>lmin.
106!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
107      DO ig=1,ngrid
108         k = lmax(ig)
109         IF (k == lmin(ig)) THEN
110            linter(ig) = k
111         ELSE
112            linter(ig) = k - zw2(ig,k) / (zw2(ig,k+1) - zw2(ig,k))
113         ENDIF
114         zmax(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig))               &
115         &        * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig)))
116      ENDDO
117     
118!===============================================================================
119! Thermal plume maximal speed and inversion layer
120!===============================================================================
121     
122!-------------------------------------------------------------------------------
123! Calcul de lmix et wmax
124!-------------------------------------------------------------------------------
125     
126      DO k=1,nlay
127         DO ig=1,ngrid
128            IF((k <= lmax(ig)).and.(k > lmin(ig))) THEN
129               IF (zw2(ig,k) < 0.) THEN
130                  print *, 'ERROR: zw2 has negative value(s)!'
131                  print *, 'ig,k', ig, k
132                  print *, 'zw2', zw2(ig,k)
133                  CALL abort
134               ENDIF
135               
136!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
137! AB: WARNING zw2 becomes its square root!
138!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
139               zw2(ig,k) = sqrt(zw2(ig,k))
140               
141               IF (zw2(ig,k) > wmax(ig)) THEN
142                  wmax(ig)  = zw2(ig,k)
143                  lmix(ig)  = k
144               ENDIF
145            ELSE
146               zw2(ig,k) = 0.
147            ENDIF
148         ENDDO
149      ENDDO
150     
151!-------------------------------------------------------------------------------
152! Calcul de zmix continu (profil parabolique des vitesses)
153!-------------------------------------------------------------------------------
154     
155      DO ig=1,ngrid
156         IF (lmix(ig) > lmin(ig)) THEN
157            IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))                        &
158            &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))             &
159            &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                   &
160            &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))) > 1e-10)   &
161            &        THEN
162               zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))              &
163               &        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)    &
164               &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))                &
165               &        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))   &
166               &        /(2.*((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))))))
170            ELSE
171               zmix(ig) = zlev(ig,lmix(ig))
172               print *, 'WARNING: problematic zmix value!'
173            ENDIF
174         ELSE
175            zmix(ig) = 0.
176         ENDIF
177      ENDDO
178     
179!-------------------------------------------------------------------------------
180! Check consistency between zmax and zmix
181!-------------------------------------------------------------------------------
182     
183      DO ig=1,ngrid
184         IF ((zmax(ig)-zmix(ig)) < 0.) THEN
185            print *, 'WARNING: we have zmix > zmax!'
186            print *, 'zmax,zmix_old,zmix_new', zmax(ig), zmix(ig), 0.9 * zmax(ig)
187            zmix(ig) = 0.9 * zmax(ig)
188            DO k=1,nlay
189               IF ((zmix(ig) >= zlev(ig,k)).and.(zmix(ig) < zlev(ig,k+1))) THEN
190                  lmix(ig) = k
191               ENDIF
192            ENDDO
193         ENDIF
194      ENDDO
195     
196     
197RETURN
198END
Note: See TracBrowser for help on using the repository browser.