source: LMDZ6/branches/blowing_snow/libf/phylmd/thermcell_height.F90 @ 5021

Last change on this file since 5021 was 4094, checked in by fhourdin, 3 years ago

Nettoyage thermiques (suite) et replay1d

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.7 KB
Line 
1      SUBROUTINE thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,  &
2     &           zw2,zlev,lmax,zmax,zmax0,zmix,wmax)
3      IMPLICIT NONE
4
5!-----------------------------------------------------------------------------
6!thermcell_height: calcul des caracteristiques du thermique: zmax,wmax,zmix
7!-----------------------------------------------------------------------------
8
9! arguments
10
11! Entree
12      integer, intent(in) :: ngrid,nlay
13      real, intent(in), dimension(ngrid) :: linter
14      real, intent(in), dimension(ngrid,nlay+1) :: zlev
15! Sortie
16      real, intent(out), dimension(ngrid) :: wmax,zmax,zmax0,zmix
17      integer, intent(out), dimension(ngrid) :: lmax
18! Les deux
19     integer, intent(inout), dimension(ngrid) :: lmix,lalim,lmin
20     real, intent(inout), dimension(ngrid,nlay+1) :: zw2
21
22! local
23     real, dimension(ngrid) :: num,denom,zlevinter
24     integer ig,l
25
26!calcul de la hauteur max du thermique
27      do ig=1,ngrid
28         lmax(ig)=lalim(ig)
29      enddo
30      do ig=1,ngrid
31         do l=nlay,lalim(ig)+1,-1
32            if (zw2(ig,l).le.1.e-10) then
33               lmax(ig)=l-1
34            endif
35         enddo
36      enddo
37
38! On traite le cas particulier qu'il faudrait éviter ou le thermique
39! atteind le haut du modele ...
40      do ig=1,ngrid
41      if ( zw2(ig,nlay) > 1.e-10 ) then
42          print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
43          lmax(ig)=nlay
44      endif
45      enddo
46
47! pas de thermique si couche 1 stable
48      do ig=1,ngrid
49         if (lmin(ig).gt.1) then
50             lmax(ig)=1
51             lmin(ig)=1
52             lalim(ig)=1
53         endif
54      enddo
55!   
56! Determination de zw2 max
57      do ig=1,ngrid
58         wmax(ig)=0.
59      enddo
60
61      do l=1,nlay
62         do ig=1,ngrid
63            if (l.le.lmax(ig)) then
64                if (zw2(ig,l).lt.0.)then
65                  print*,'pb2 zw2<0'
66                endif
67                zw2(ig,l)=sqrt(zw2(ig,l))
68                wmax(ig)=max(wmax(ig),zw2(ig,l))
69            else
70                 zw2(ig,l)=0.
71            endif
72          enddo
73      enddo
74
75!   Longueur caracteristique correspondant a la hauteur des thermiques.
76      do  ig=1,ngrid
77         zmax(ig)=0.
78         zlevinter(ig)=zlev(ig,1)
79      enddo
80
81!     if (iflag_thermals_ed.ge.1) then
82      if (1==0) then
83!CR:date de quand le calcul du zmax continu etait buggue
84         num(:)=0.
85         denom(:)=0.
86         do ig=1,ngrid
87          do l=1,nlay
88             num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
89             denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
90          enddo
91       enddo
92       do ig=1,ngrid
93       if (denom(ig).gt.1.e-10) then
94          zmax(ig)=2.*num(ig)/denom(ig)
95          zmax0(ig)=zmax(ig)
96       endif
97       enddo
98 
99      else
100!CR:Calcul de zmax continu via le linter     
101      do  ig=1,ngrid
102! calcul de zlevinter
103          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
104     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
105     &    -zlev(ig,lmax(ig)))
106!pour le cas ou on prend tjs lmin=1
107!       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
108       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
109       zmax0(ig)=zmax(ig)
110      enddo
111
112
113      endif
114!endif iflag_thermals_ed
115!
116! def de  zmix continu (profil parabolique des vitesses)
117      do ig=1,ngrid
118           if (lmix(ig).gt.1) then
119! test
120              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
121     &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))  &
122     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
123     &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)  &
124     &        then
125!             
126            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
127     &        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)  &
128     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
129     &        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))  &
130     &        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
131     &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))  &
132     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
133     &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
134              else
135              zmix(ig)=zlev(ig,lmix(ig))
136              print*,'pb zmix'
137              endif
138          else
139              zmix(ig)=0.
140          endif
141!test
142         if ((zmax(ig)-zmix(ig)).le.0.) then
143            zmix(ig)=0.9*zmax(ig)
144!            print*,'pb zmix>zmax'
145         endif
146      enddo
147!
148! calcul du nouveau lmix correspondant
149      do ig=1,ngrid
150         do l=1,nlay
151            if (zmix(ig).ge.zlev(ig,l).and.  &
152     &          zmix(ig).lt.zlev(ig,l+1)) then
153              lmix(ig)=l
154             endif
155          enddo
156      enddo
157!
158 RETURN
159      end
Note: See TracBrowser for help on using the repository browser.