GCC Code Coverage Report


Directory: ./
File: phys/thermcell_dv2.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 70 0.0%
Branches: 0 56 0.0%

Line Branch Exec Source
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 !----------------------
103 do 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 !--------------------
137 enddo
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
194