source: LMDZ5/trunk/libf/phylmd/thermcell_height.F90 @ 4404

Last change on this file since 4404 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

  • moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
  • created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
  • created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
  • reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

EM
  • 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,lev_out)                           
3
4!-----------------------------------------------------------------------------
5!thermcell_height: calcul des caracteristiques du thermique: zmax,wmax,zmix
6!-----------------------------------------------------------------------------
7      IMPLICIT NONE
8#include "thermcell.h"
9
10      INTEGER ig,l
11      INTEGER ngrid,nlay
12      INTEGER lalim(ngrid),lmin(ngrid)
13      INTEGER lmix(ngrid)
14      REAL linter(ngrid)
15      integer lev_out                           ! niveau pour les print
16
17      REAL zw2(ngrid,nlay+1)
18      REAL zlev(ngrid,nlay+1)
19
20      REAL wmax(ngrid)
21      INTEGER lmax(ngrid)
22      REAL zmax(ngrid)
23      REAL zmax0(ngrid)
24      REAL zmix(ngrid)
25      REAL num(ngrid)
26      REAL denom(ngrid)
27
28      REAL zlevinter(ngrid)
29
30!calcul de la hauteur max du thermique
31      do ig=1,ngrid
32         lmax(ig)=lalim(ig)
33      enddo
34      do ig=1,ngrid
35         do l=nlay,lalim(ig)+1,-1
36            if (zw2(ig,l).le.1.e-10) then
37               lmax(ig)=l-1
38            endif
39         enddo
40      enddo
41
42! On traite le cas particulier qu'il faudrait éviter ou le thermique
43! atteind le haut du modele ...
44      do ig=1,ngrid
45      if ( zw2(ig,nlay) > 1.e-10 ) then
46          print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
47          lmax(ig)=nlay
48      endif
49      enddo
50
51! pas de thermique si couche 1 stable
52      do ig=1,ngrid
53         if (lmin(ig).gt.1) then
54             lmax(ig)=1
55             lmin(ig)=1
56             lalim(ig)=1
57         endif
58      enddo
59!   
60! Determination de zw2 max
61      do ig=1,ngrid
62         wmax(ig)=0.
63      enddo
64
65      do l=1,nlay
66         do ig=1,ngrid
67            if (l.le.lmax(ig)) then
68                if (zw2(ig,l).lt.0.)then
69                  print*,'pb2 zw2<0'
70                endif
71                zw2(ig,l)=sqrt(zw2(ig,l))
72                wmax(ig)=max(wmax(ig),zw2(ig,l))
73            else
74                 zw2(ig,l)=0.
75            endif
76          enddo
77      enddo
78
79!   Longueur caracteristique correspondant a la hauteur des thermiques.
80      do  ig=1,ngrid
81         zmax(ig)=0.
82         zlevinter(ig)=zlev(ig,1)
83      enddo
84
85!     if (iflag_thermals_ed.ge.1) then
86      if (1==0) then
87!CR:date de quand le calcul du zmax continu etait buggue
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!CR:Calcul de zmax continu via le linter     
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.