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

Last change on this file since 3058 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
RevLine 
[2060]1!
2!
3!
[2480]4SUBROUTINE thermcell_height(ngrid,nlay,lmin,lmix,lmax,                        &
5                            zlev,zmin,zmix,zmax,zw2,wmax,f_star)
[2060]6     
7     
[2127]8!===============================================================================
9!  Purpose: calcul des caracteristiques du thermique: zmax,wmax,zmix
10!===============================================================================
[2060]11     
12      IMPLICIT NONE
13     
14     
[2127]15!===============================================================================
[2060]16! Declaration
[2127]17!===============================================================================
[2060]18     
[2127]19!     Inputs:
20!     -------
[2060]21     
[2177]22      INTEGER, INTENT(in) :: ngrid
23      INTEGER, INTENT(in) :: nlay
24      INTEGER, INTENT(in) :: lmin(ngrid)
[2060]25     
[2177]26      REAL, INTENT(in) :: zlev(ngrid,nlay+1)
27      REAL, INTENT(in) :: f_star(ngrid,nlay+1)
[2060]28     
[2127]29!     Outputs:
30!     --------
[2060]31     
[2177]32      INTEGER, INTENT(out) :: lmix(ngrid)
33      INTEGER, INTENT(out) :: lmax(ngrid)
[2060]34     
[2177]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
[2060]39     
[2177]40      REAL, INTENT(inout) :: zw2(ngrid,nlay+1)        ! Square vertical speed (becomes its square root)
41     
[2127]42!     Local:
43!     ------
[2060]44     
[2177]45      INTEGER ig, k
[2060]46     
[2177]47      REAL linter(ngrid)                              ! Level (continuous) of maximal vertical speed
[2060]48     
[2127]49!===============================================================================
[2060]50! Initialization
[2127]51!===============================================================================
[2060]52     
[2177]53      linter(:) = lmin(:)
54      lmix(:) = lmin(:)
55      lmax(:) = lmin(:)
[2060]56     
[2177]57      wmax(:) = 0.
[2060]58     
[2177]59      zmin(:) = 0.
60      zmix(:) = 0.
61      zmax(:) = 0.
[2060]62     
[2127]63!===============================================================================
[2143]64! Thermal plume height
[2127]65!===============================================================================
[2060]66     
[2127]67!-------------------------------------------------------------------------------
[2143]68! Calcul de zmin
69!-------------------------------------------------------------------------------
70     
[2177]71      DO ig=1,ngrid
72         zmin(ig) = zlev(ig,lmin(ig))
[2143]73      ENDDO
74     
75!-------------------------------------------------------------------------------
[2060]76! Calcul de lmax
[2127]77!-------------------------------------------------------------------------------
[2060]78     
79      DO ig=1,ngrid
[2177]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
[2060]83            ENDIF
84         ENDDO
85      ENDDO
86     
[2127]87!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2143]88! AB: Problematic case where thermal plume reaches the top of the model...
[2127]89!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]90      DO ig=1,ngrid
[2143]91         IF (zw2(ig,nlay) > 0.) THEN
[2127]92            print *, 'WARNING: a thermal plume reaches the highest layer!'
[2177]93            print *, 'ig,k', ig, nlay
[2127]94            print *, 'zw2', zw2(ig,nlay)
95            lmax(ig) = nlay
[2060]96         ENDIF
97      ENDDO
98     
[2127]99!-------------------------------------------------------------------------------
[2177]100! Calcul de zmax continu via linter
[2127]101!-------------------------------------------------------------------------------
[2060]102     
[2127]103!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2177]104! AB: lmin=lmax means the plume is not active and then zw2=0 at each level.
105!     Otherwise we have lmax>lmin.
[2127]106!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]107      DO ig=1,ngrid
[2177]108         k = lmax(ig)
109         IF (k == lmin(ig)) THEN
110            linter(ig) = k
[2060]111         ELSE
[2177]112            linter(ig) = k - zw2(ig,k) / (zw2(ig,k+1) - zw2(ig,k))
[2060]113         ENDIF
[2177]114         zmax(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig))               &
115         &        * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig)))
[2060]116      ENDDO
117     
[2127]118!===============================================================================
[2143]119! Thermal plume maximal speed and inversion layer
[2127]120!===============================================================================
[2060]121     
[2127]122!-------------------------------------------------------------------------------
[2060]123! Calcul de lmix et wmax
[2127]124!-------------------------------------------------------------------------------
[2060]125     
[2177]126      DO k=1,nlay
[2060]127         DO ig=1,ngrid
[2177]128            IF((k <= lmax(ig)).and.(k > lmin(ig))) THEN
129               IF (zw2(ig,k) < 0.) THEN
[2060]130                  print *, 'ERROR: zw2 has negative value(s)!'
[2177]131                  print *, 'ig,k', ig, k
132                  print *, 'zw2', zw2(ig,k)
[2060]133                  CALL abort
134               ENDIF
135               
[2127]136!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2177]137! AB: WARNING zw2 becomes its square root!
[2127]138!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2177]139               zw2(ig,k) = sqrt(zw2(ig,k))
[2060]140               
[2177]141               IF (zw2(ig,k) > wmax(ig)) THEN
142                  wmax(ig)  = zw2(ig,k)
143                  lmix(ig)  = k
[2060]144               ENDIF
145            ELSE
[2177]146               zw2(ig,k) = 0.
[2060]147            ENDIF
148         ENDDO
149      ENDDO
150     
[2127]151!-------------------------------------------------------------------------------
[2060]152! Calcul de zmix continu (profil parabolique des vitesses)
[2127]153!-------------------------------------------------------------------------------
[2060]154     
155      DO ig=1,ngrid
[2143]156         IF (lmix(ig) > lmin(ig)) THEN
[2060]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))                   &
[2143]160            &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))) > 1e-10)   &
[2060]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))                &
[2143]169               &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
[2060]170            ELSE
[2143]171               zmix(ig) = zlev(ig,lmix(ig))
[2127]172               print *, 'WARNING: problematic zmix value!'
[2060]173            ENDIF
174         ELSE
175            zmix(ig) = 0.
176         ENDIF
177      ENDDO
178     
[2143]179!-------------------------------------------------------------------------------
[2177]180! Check consistency between zmax and zmix
[2143]181!-------------------------------------------------------------------------------
182     
[2060]183      DO ig=1,ngrid
[2143]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)
[2060]187            zmix(ig) = 0.9 * zmax(ig)
[2177]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
[2060]193         ENDIF
194      ENDDO
195     
196     
197RETURN
198END
Note: See TracBrowser for help on using the repository browser.