source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/thermcell_dv2.F90 @ 3158

Last change on this file since 3158 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: 5.7 KB
Line 
1      subroutine thermcell_dv2(ngrid,nlay,ptimestep,fm,entr,masse  &
2     &    ,fraca,larga  &
3     &    ,u,v,du,dv,ua,va,lev_out)
4      USE print_control_mod, ONLY: prt_level,lunout
5      implicit none
6
7!=======================================================================
8!
9!   Calcul du transport verticale dans la couche limite en presence
10!   de "thermiques" explicitement representes
11!   calcul du dq/dt une fois qu'on connait les ascendances
12!
13! Vectorisation, FH : 2010/03/08
14!
15!=======================================================================
16
17
18      integer ngrid,nlay
19
20      real ptimestep
21      real masse(ngrid,nlay),fm(ngrid,nlay+1)
22      real fraca(ngrid,nlay+1)
23      real larga(ngrid)
24      real entr(ngrid,nlay)
25      real u(ngrid,nlay)
26      real ua(ngrid,nlay)
27      real du(ngrid,nlay)
28      real v(ngrid,nlay)
29      real va(ngrid,nlay)
30      real dv(ngrid,nlay)
31      integer lev_out                           ! niveau pour les print
32
33      real qa(ngrid,nlay),detr(ngrid,nlay),zf,zf2
34      real wvd(ngrid,nlay+1),wud(ngrid,nlay+1)
35      real gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1)
36      real ue(ngrid,nlay),ve(ngrid,nlay)
37      LOGICAL ltherm(ngrid,nlay)
38      real dua(ngrid,nlay),dva(ngrid,nlay)
39      integer iter
40
41      integer ig,k,nlarga0
42
43!-------------------------------------------------------------------------
44
45!   calcul du detrainement
46!---------------------------
47
48!      print*,'THERMCELL DV2 OPTIMISE 3'
49
50      nlarga0=0.
51
52      do k=1,nlay
53         do ig=1,ngrid
54            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
55         enddo
56      enddo
57
58!   calcul de la valeur dans les ascendances
59      do ig=1,ngrid
60         ua(ig,1)=u(ig,1)
61         va(ig,1)=v(ig,1)
62         ue(ig,1)=u(ig,1)
63         ve(ig,1)=v(ig,1)
64      enddo
65
66      IF(prt_level>9)WRITE(lunout,*)                                    &
67     &      'WARNING on initialise gamma(1:ngrid,1)=0.'
68      gamma(1:ngrid,1)=0.
69      do k=2,nlay
70         do ig=1,ngrid
71            ltherm(ig,k)=(fm(ig,k+1)+detr(ig,k))*ptimestep > 1.e-5*masse(ig,k)
72            if(ltherm(ig,k).and.larga(ig)>0.) then
73               gamma0(ig,k)=masse(ig,k)  &
74     &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
75     &         *0.5/larga(ig)  &
76     &         *1.
77            else
78               gamma0(ig,k)=0.
79            endif
80            if (ltherm(ig,k).and.larga(ig)<=0.) nlarga0=nlarga0+1
81         enddo
82      enddo
83
84      gamma(:,:)=0.
85
86      do k=2,nlay
87
88         do ig=1,ngrid
89            if (ltherm(ig,k)) then
90               dua(ig,k)=ua(ig,k-1)-u(ig,k-1)
91               dva(ig,k)=va(ig,k-1)-v(ig,k-1)
92            else
93               ua(ig,k)=u(ig,k)
94               va(ig,k)=v(ig,k)
95               ue(ig,k)=u(ig,k)
96               ve(ig,k)=v(ig,k)
97            endif
98         enddo
99
100
101! Debut des iterations
102!----------------------
103do iter=1,5
104         do ig=1,ngrid
105! Pour memoire : calcul prenant en compte la fraction reelle
106!              zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
107!              zf2=1./(1.-zf)
108! Calcul avec fraction infiniement petite
109               zf=0.
110               zf2=1.
111
112!  la première fois on multiplie le coefficient de freinage
113!  par le module du vent dans la couche en dessous.
114!  Mais pourquoi donc ???
115               if (ltherm(ig,k)) then
116!   On choisit une relaxation lineaire.
117!                 gamma(ig,k)=gamma0(ig,k)
118!   On choisit une relaxation quadratique.
119                  gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2)
120                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)  &
121     &               +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))  &
122     &               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2  &
123     &                 +gamma(ig,k))
124                  va(ig,k)=(fm(ig,k)*va(ig,k-1)  &
125     &               +(zf2*entr(ig,k)+gamma(ig,k))*v(ig,k))  &
126     &               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2  &
127     &                 +gamma(ig,k))
128!                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua(ig,k),dva(ig,k)
129                  dua(ig,k)=ua(ig,k)-u(ig,k)
130                  dva(ig,k)=va(ig,k)-v(ig,k)
131                  ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
132                  ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
133               endif
134         enddo
135! Fin des iterations
136!--------------------
137enddo
138
139      enddo ! k=2,nlay
140
141
142! Calcul du flux vertical de moment dans l'environnement.
143!---------------------------------------------------------
144      do k=2,nlay
145         do ig=1,ngrid
146            wud(ig,k)=fm(ig,k)*ue(ig,k)
147            wvd(ig,k)=fm(ig,k)*ve(ig,k)
148         enddo
149      enddo
150      do ig=1,ngrid
151         wud(ig,1)=0.
152         wud(ig,nlay+1)=0.
153         wvd(ig,1)=0.
154         wvd(ig,nlay+1)=0.
155      enddo
156
157! calcul des tendances.
158!-----------------------
159      do k=1,nlay
160         do ig=1,ngrid
161            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)  &
162     &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
163     &               -wud(ig,k)+wud(ig,k+1))  &
164     &               /masse(ig,k)
165            dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)  &
166     &               -(entr(ig,k)+gamma(ig,k))*ve(ig,k)  &
167     &               -wvd(ig,k)+wvd(ig,k+1))  &
168     &               /masse(ig,k)
169         enddo
170      enddo
171
172
173! Sorties eventuelles.
174!----------------------
175
176   if(prt_level.GE.10) then
177      do k=1,nlay
178         do ig=1,ngrid
179           print*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), &
180     &   entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), &
181     &   masse(ig,k)
182         enddo
183      enddo
184   endif
185!
186     if (nlarga0>0) then
187          print*,'WARNING !!!!!! DANS THERMCELL_DV2 '
188          print*,nlarga0,' points pour lesquels laraga=0. dans un thermique'
189          print*,'Il faudrait decortiquer ces points'
190     endif
191
192      return
193      end
Note: See TracBrowser for help on using the repository browser.