source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/thermcell_height.F90 @ 4104

Last change on this file since 4104 was 1026, checked in by lmdzadmin, 16 years ago

Modifs thermiques
FH/IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.3 KB
Line 
1      SUBROUTINE thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,  &
2     &           zw2,zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)                           
3
4!-----------------------------------------------------------------------------
5!thermcell_height: calcul des caracteristiques du thermique: zmax,wmax,zmix
6!-----------------------------------------------------------------------------
7      IMPLICIT NONE
8#include "iniprint.h"
9#include "thermcell.h"
10
11      INTEGER ig,l
12      INTEGER ngrid,nlay
13      INTEGER lalim(ngrid),lmin(ngrid)
14      INTEGER lmix(ngrid)
15      REAL linter(ngrid)
16      integer lev_out                           ! niveau pour les print
17
18      REAL zw2(ngrid,nlay+1)
19      REAL zlev(ngrid,nlay+1)
20
21      REAL wmax(ngrid)
22      INTEGER lmax(ngrid)
23      REAL zmax(ngrid)
24      REAL zmax0(ngrid)
25      REAL zmix(ngrid)
26      REAL num(ngrid)
27      REAL denom(ngrid)
28
29      REAL zlevinter(ngrid)
30
31!calcul de la hauteur max du thermique
32      do ig=1,ngrid
33         lmax(ig)=lalim(ig)
34      enddo
35      do ig=1,ngrid
36         do l=nlay,lalim(ig)+1,-1
37            if (zw2(ig,l).le.1.e-10) then
38               lmax(ig)=l-1
39            endif
40         enddo
41      enddo
42! pas de thermique si couche 1 stable
43      do ig=1,ngrid
44         if (lmin(ig).gt.1) then
45             lmax(ig)=1
46             lmin(ig)=1
47             lalim(ig)=1
48         endif
49      enddo
50!   
51! Determination de zw2 max
52      do ig=1,ngrid
53         wmax(ig)=0.
54      enddo
55
56      do l=1,nlay
57         do ig=1,ngrid
58            if (l.le.lmax(ig)) then
59                if (zw2(ig,l).lt.0.)then
60                  print*,'pb2 zw2<0'
61                endif
62                zw2(ig,l)=sqrt(zw2(ig,l))
63                wmax(ig)=max(wmax(ig),zw2(ig,l))
64            else
65                 zw2(ig,l)=0.
66            endif
67          enddo
68      enddo
69
70!   Longueur caracteristique correspondant a la hauteur des thermiques.
71      do  ig=1,ngrid
72         zmax(ig)=0.
73         zlevinter(ig)=zlev(ig,1)
74      enddo
75
76      if (iflag_thermals_ed.ge.1) then
77
78         num(:)=0.
79         denom(:)=0.
80         do ig=1,ngrid
81          do l=1,nlay
82             num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
83             denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
84          enddo
85       enddo
86       do ig=1,ngrid
87       if (denom(ig).gt.1.e-10) then
88          zmax(ig)=2.*num(ig)/denom(ig)
89          zmax0(ig)=zmax(ig)
90       endif
91       enddo
92
93       else
94
95      do  ig=1,ngrid
96! calcul de zlevinter
97          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
98     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
99     &    -zlev(ig,lmax(ig)))
100!pour le cas ou on prend tjs lmin=1
101!       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
102       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
103       zmax0(ig)=zmax(ig)
104      enddo
105
106
107      endif
108!endif iflag_thermals_ed
109!
110! def de  zmix continu (profil parabolique des vitesses)
111      do ig=1,ngrid
112           if (lmix(ig).gt.1) then
113! test
114              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
115     &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))  &
116     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
117     &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)  &
118     &        then
119!             
120            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
121     &        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)  &
122     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
123     &        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))  &
124     &        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
125     &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))  &
126     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
127     &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
128              else
129              zmix(ig)=zlev(ig,lmix(ig))
130              print*,'pb zmix'
131              endif
132          else
133              zmix(ig)=0.
134          endif
135!test
136         if ((zmax(ig)-zmix(ig)).le.0.) then
137            zmix(ig)=0.9*zmax(ig)
138!            print*,'pb zmix>zmax'
139         endif
140      enddo
141!
142! calcul du nouveau lmix correspondant
143      do ig=1,ngrid
144         do l=1,nlay
145            if (zmix(ig).ge.zlev(ig,l).and.  &
146     &          zmix(ig).lt.zlev(ig,l+1)) then
147              lmix(ig)=l
148             endif
149          enddo
150      enddo
151!
152      return
153      end
Note: See TracBrowser for help on using the repository browser.