source: LMDZ5/trunk/libf/phylmd/cv3_vertmix.F @ 1912

Last change on this file since 1912 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • 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: 5.0 KB
Line 
1      SUBROUTINE cv3_vertmix(len,nd,iflag,plim1,plim2,p,ph,t,q,u,v
2     :                     ,w,wi,nk,tmix,thmix,qmix,qsmix
3     :                     ,umix,vmix,plcl)
4***************************************************************
5*                                                             *
6* CV3_VERTMIX   Brassage adiabatique d'une couche d'epaisseur *
7*               arbitraire.                                   *
8*                                                             *
9* written by   : Grandpeix Jean-Yves, 28/12/2001, 13.14.24    *
10* modified by :  Filiberti M-A 06/2005 vectorisation          *
11***************************************************************
12*
13       implicit none
14C==============================================================
15C
16C vertmix : determine theta et r du melange obtenu en brassant
17C adiabatiquement entre plim1 et plim2, avec une ponderation w.
18C
19C===============================================================
20
21      include "cvthermo.h"
22      include "YOETHF.h"
23      include "YOMCST.h"
24      include "FCTTRE.h"
25c input :
26      integer nd,len
27      integer nk(len),iflag(len)
28      real t(len,nd),q(len,nd),w(nd)
29      real u(len,nd),v(len,nd)
30      real p(len,nd),ph(len,nd+1)
31      real plim1(len),plim2(len)
32c output :
33      real tmix(len),thmix(len),qmix(len),wi(len,nd)
34      real umix(len),vmix(len)
35      real qsmix(len)
36      real plcl(len)
37c internal variables :
38      integer j1(len),j2(len),niflag7
39      real A,B
40      real ahm(len),dpw(len),coef(len)
41      real p1(len,nd),p2(len,nd)
42      real rdcp(len),a2(len),b2(len),pnk(len)
43      real rh(len),chi(len)
44      real cpn
45      real x,y,p0,p0m1,zdelta,zcor
46
47      integer i,j
48
49      do j = 1,nd
50        do i=1,len
51          if (plim1(i).le.ph(i,j)) j1(i) = j
52          if (plim2(i).ge.ph(i,j+1).and.plim2(i).lt.ph(i,j)) j2(i) = j
53        enddo
54      enddo
55c
56      do j=1,nd
57        do i = 1,len
58          wi(i,j) = 0.
59        enddo
60      enddo
61      do i = 1,len
62       ahm(i)=0.
63       qmix(i)=0.
64       umix(i)=0.
65       vmix(i)=0.
66       dpw(i) =0.
67       a2(i)=0.0
68       b2(i) = 0.
69       pnk(i) = p(i,nk(i))
70      enddo
71c
72      p0 = 1000.
73      p0m1 = 1./p0
74c
75      do i=1,len
76        coef(i) = 1./(plim1(i)-plim2(i))
77      end do
78c
79      do  j=1,nd
80        do i=1,len
81          if (j.ge.j1(i).and.j.le.j2(i)) then
82            p1(i,j) = min(ph(i,j),plim1(i))
83            p2(i,j) = max(ph(i,j+1),plim2(i))
84cCRtest:couplage thermiques: deja normalise
85c             wi(i,j) = w(j)
86c             print*,'wi',wi(i,j)
87            wi(i,j) = w(j)*(p1(i,j)-p2(i,j))*coef(i)
88            dpw(i) = dpw(i)+wi(i,j)
89          endif
90        end do
91      end do
92cCR:print
93c      do i=1,len
94c         print*,'plim',plim1(i),plim2(i)
95c      enddo
96      do  j=1,nd
97        do i=1,len
98          if (j.ge.j1(i).and.j.le.j2(i)) then
99            wi(i,j)=wi(i,j)/dpw(i)
100            ahm(i)=ahm(i)+(cpd*(1.-q(i,j))+q(i,j)*cpv)*t(i,j)*wi(i,j)
101            qmix(i)=qmix(i)+q(i,j)*wi(i,j)
102            umix(i)=umix(i)+u(i,j)*wi(i,j)
103            vmix(i)=vmix(i)+v(i,j)*wi(i,j)
104          endif
105        end do
106      end do
107c
108      do i=1,len
109         rdcp(i)=(rrd*(1.-qmix(i))+qmix(i)*rrv)/
110     :            (cpd*(1.-qmix(i))+qmix(i)*cpv)
111      end do
112c
113
114c
115      do 20 j=1,nd
116        do 18 i=1,len
117          if (j.ge.j1(i).and.j.le.j2(i)) then
118cc            x=(.5*(p1(i,j)+p2(i,j))*p0m1)**rdcp(i)
119            y=(.5*(p1(i,j)+p2(i,j))/pnk(i))**rdcp(i)
120cc            a2(i)=a2(i)+(cpd*(1.-qmix(i))+qmix(i)*cpv)*x*wi(i,j)
121            b2(i)=b2(i)+(cpd*(1.-qmix(i))+qmix(i)*cpv)*y*wi(i,j)
122          endif
123   18  continue
124   20  continue
125c
126       do i=1,len
127         tmix(i) = ahm(i)/b2(i)
128         thmix(i) =tmix(i)*(p0/pnk(i))**rdcp(i)
129c         print*,'thmix ahm',ahm(i),b2(i)
130c         print*,'thmix t',tmix(i),p0
131c         print*,'thmix p',pnk(i),rdcp(i)
132c         print*,'thmix',thmix(i)
133cc         thmix(i) = ahm(i)/a2(i)
134cc         tmix(i)= thmix(i)*(pnk(i)*p0m1)**rdcp(i)
135         zdelta=max(0.,sign(1.,rtt-tmix(i)))
136         qsmix(i)= r2es*FOEEW(tmix(i),zdelta)/(pnk(i)*100.)
137         qsmix(i)=min(0.5,qsmix(i))
138         zcor=1./(1.-retv*qsmix(i))
139         qsmix(i)=qsmix(i)*zcor
140       end do
141c
142!-------------------------------------------------------------------
143! --- Calculate lifted condensation level of air at parcel origin level
144! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
145!-------------------------------------------------------------------
146
147       A = 1669.0 ! convect3
148       B = 122.0  ! convect3
149
150
151       niflag7=0
152       do 260 i=1,len
153
154        if (iflag(i).ne.7) then ! modif sb Jun7th 2002
155c
156        rh(i)=qmix(i)/qsmix(i)
157        chi(i)=tmix(i)/(A-B*rh(i)-tmix(i)) ! convect3
158c   ATTENTION, la LIGNE DESSOUS A ETE RAJOUTEE ARBITRAIREMENT ET
159c  MASQUE UN PB POTENTIEL
160        chi(i)=max(chi(i),0.)
161        rh(i)=max(rh(i),0.)
162        plcl(i)=pnk(i)*(rh(i)**chi(i))
163        if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
164     &   .and.(iflag(i).eq.0))iflag(i)=8
165
166        else
167
168          niflag7=niflag7+1
169          plcl(i)=plim2(i)
170c
171        endif ! iflag=7
172
173c      print*,'NIFLAG7  =',niflag7
174
175 260   continue
176
177      return
178      end
179
Note: See TracBrowser for help on using the repository browser.