source: LMDZ4/branches/LMDZ4V5.0-LF/libf/phylmd/thermcell_height.F90 @ 3603

Last change on this file since 3603 was 1311, checked in by Laurent Fairhead, 15 years ago

Modifications to thermals for TKE transport

Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.6 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
43! On traite le cas particulier qu'il faudrait éviter ou le thermique
44! atteind le haut du modele ...
45      do ig=1,ngrid
46      if ( zw2(ig,nlay) > 1.e-10 ) then
47          print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
48          lmax(ig)=nlay
49      endif
50      enddo
51
52! pas de thermique si couche 1 stable
53      do ig=1,ngrid
54         if (lmin(ig).gt.1) then
55             lmax(ig)=1
56             lmin(ig)=1
57             lalim(ig)=1
58         endif
59      enddo
60!   
61! Determination de zw2 max
62      do ig=1,ngrid
63         wmax(ig)=0.
64      enddo
65
66      do l=1,nlay
67         do ig=1,ngrid
68            if (l.le.lmax(ig)) then
69                if (zw2(ig,l).lt.0.)then
70                  print*,'pb2 zw2<0'
71                endif
72                zw2(ig,l)=sqrt(zw2(ig,l))
73                wmax(ig)=max(wmax(ig),zw2(ig,l))
74            else
75                 zw2(ig,l)=0.
76            endif
77          enddo
78      enddo
79
80!   Longueur caracteristique correspondant a la hauteur des thermiques.
81      do  ig=1,ngrid
82         zmax(ig)=0.
83         zlevinter(ig)=zlev(ig,1)
84      enddo
85
86      if (iflag_thermals_ed.ge.1) then
87
88         num(:)=0.
89         denom(:)=0.
90         do ig=1,ngrid
91          do l=1,nlay
92             num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
93             denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
94          enddo
95       enddo
96       do ig=1,ngrid
97       if (denom(ig).gt.1.e-10) then
98          zmax(ig)=2.*num(ig)/denom(ig)
99          zmax0(ig)=zmax(ig)
100       endif
101       enddo
102
103       else
104
105      do  ig=1,ngrid
106! calcul de zlevinter
107          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
108     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
109     &    -zlev(ig,lmax(ig)))
110!pour le cas ou on prend tjs lmin=1
111!       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
112       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
113       zmax0(ig)=zmax(ig)
114      enddo
115
116
117      endif
118!endif iflag_thermals_ed
119!
120! def de  zmix continu (profil parabolique des vitesses)
121      do ig=1,ngrid
122           if (lmix(ig).gt.1) then
123! test
124              if (((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))))).gt.1e-10)  &
128     &        then
129!             
130            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
131     &        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)  &
132     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
133     &        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))  &
134     &        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
135     &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))  &
136     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
137     &        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
138              else
139              zmix(ig)=zlev(ig,lmix(ig))
140              print*,'pb zmix'
141              endif
142          else
143              zmix(ig)=0.
144          endif
145!test
146         if ((zmax(ig)-zmix(ig)).le.0.) then
147            zmix(ig)=0.9*zmax(ig)
148!            print*,'pb zmix>zmax'
149         endif
150      enddo
151!
152! calcul du nouveau lmix correspondant
153      do ig=1,ngrid
154         do l=1,nlay
155            if (zmix(ig).ge.zlev(ig,l).and.  &
156     &          zmix(ig).lt.zlev(ig,l+1)) then
157              lmix(ig)=l
158             endif
159          enddo
160      enddo
161!
162      return
163      end
Note: See TracBrowser for help on using the repository browser.