/[lmdze]/trunk/phylmd/clqh.f
ViewVC logotype

Contents of /trunk/phylmd/clqh.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 202 - (show annotations)
Wed Jun 8 12:23:41 2016 UTC (7 years, 10 months ago) by guez
Original Path: trunk/Sources/phylmd/clqh.f
File size: 10020 byte(s) 
Promoted lmt_pas from local variable of physiq to variable of module
conf_gcm_m.

Removed variable run_off of module interface_surf. Was not
used. Called run_off_ter in LMDZ, but not used nor printed there
either.

Simplified logic in interfoce_lim. The way it was convoluted with
interfsurf_hq and clmain was quite a mess. Extracted reading of SST
into a separate procedure: read_sst. We do not need SST and pctsrf_new
at the same time: SST is not needed for sea-ice surface. I did not
like this programming: going through the procedure repeatedly for
different purposes and testing inside whether there was something to
do or it was already done. Reading is now only controlled by itap and
lmt_pas, instead of debut, jour, jour_lu and deja_lu. Now we do not
copy from pct_tmp to pctsrf_new every time step.

Simplified processing of pctsrf in clmain and below. It was quite
troubling: pctsrf_new was intent out in interfoce_lim but only defined
for ocean and sea-ice. Also the idea of having arrays for all
surfaces, pcsrf and pctsrf_new, in interfsurf_hq, which is called for
a particular surface, was troubling. pctsrf_new for all surfaces was
intent out in intefsurf_hq, but not defined for all surfaces at each
call. Removed argument pctsrf_new of clmain: was a duplicate of pctsrf
on output, and not used in physiq. Replaced pctsrf_new in clmain by
pctsrf_new_oce and pctsrf_new_sic, which were the only ones modified.
1 module clqh_m
2
3 IMPLICIT none
4
5 contains
6
7 SUBROUTINE clqh(dtime, jour, debut, rlat, knon, nisurf, knindex, &
8 tsoil, qsol, rmu0, rugos, rugoro, u1lay, v1lay, coef, t, q, ts, paprs, &
9 pplay, delp, radsol, albedo, snow, qsurf, precip_rain, precip_snow, &
10 fder, fluxlat, pctsrf_new_sic, agesno, d_t, d_q, d_ts, z0_new, flux_t, &
11 flux_q, dflux_s, dflux_l, fqcalving, ffonte, run_off_lic_0)
12
13 ! Author: Z. X. Li (LMD/CNRS)
14 ! Date: 1993/08/18
15 ! Objet : diffusion verticale de "q" et de "h"
16
17 USE conf_phys_m, ONLY: iflag_pbl
18 USE dimphy, ONLY: klev, klon
19 USE dimsoil, ONLY: nsoilmx
20 USE indicesol, ONLY: nbsrf
21 USE interfsurf_hq_m, ONLY: interfsurf_hq
22 USE suphec_m, ONLY: rcpd, rd, rg, rkappa
23
24 REAL, intent(in):: dtime ! intervalle du temps (s)
25 integer, intent(in):: jour ! jour de l'annee en cours
26 logical, intent(in):: debut
27 real, intent(in):: rlat(klon)
28 INTEGER, intent(in):: knon
29 integer, intent(in):: nisurf
30 integer, intent(in):: knindex(:) ! (knon)
31
32 REAL tsoil(klon, nsoilmx)
33
34 REAL, intent(inout):: qsol(klon)
35 ! column-density of water in soil, in kg m-2
36
37 real, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal
38 real rugos(klon) ! rugosite
39 REAL rugoro(klon)
40 REAL u1lay(klon) ! vitesse u de la 1ere couche (m / s)
41 REAL v1lay(klon) ! vitesse v de la 1ere couche (m / s)
42
43 REAL, intent(in):: coef(:, :) ! (knon, klev)
44 ! Le coefficient d'echange (m**2 / s) multiplie par le cisaillement
45 ! du vent (dV / dz). La premiere valeur indique la valeur de Cdrag
46 ! (sans unite).
47
48 REAL t(klon, klev) ! temperature (K)
49 REAL q(klon, klev) ! humidite specifique (kg / kg)
50 REAL, intent(in):: ts(klon) ! temperature du sol (K)
51 REAL paprs(klon, klev+1) ! pression a inter-couche (Pa)
52 REAL pplay(klon, klev) ! pression au milieu de couche (Pa)
53 REAL delp(klon, klev) ! epaisseur de couche en pression (Pa)
54 REAL radsol(klon) ! ray. net au sol (Solaire+IR) W / m2
55 REAL, intent(inout):: albedo(:) ! (knon) albedo de la surface
56 REAL, intent(inout):: snow(klon) ! hauteur de neige
57 REAL qsurf(klon) ! humidite de l'air au dessus de la surface
58
59 real, intent(in):: precip_rain(klon)
60 ! liquid water mass flux (kg / m2 / s), positive down
61
62 real, intent(in):: precip_snow(klon)
63 ! solid water mass flux (kg / m2 / s), positive down
64
65 real, intent(inout):: fder(klon)
66 real fluxlat(klon)
67 real, intent(in):: pctsrf_new_sic(:) ! (klon)
68 REAL, intent(inout):: agesno(:) ! (knon)
69 REAL d_t(klon, klev) ! incrementation de "t"
70 REAL d_q(klon, klev) ! incrementation de "q"
71 REAL, intent(out):: d_ts(:) ! (knon) incrementation de "ts"
72 real z0_new(klon)
73 REAL flux_t(klon, klev) ! (diagnostic) flux de la chaleur
74 ! sensible, flux de Cp*T, positif vers
75 ! le bas: j / (m**2 s) c.a.d.: W / m2
76 REAL flux_q(klon, klev) ! flux de la vapeur d'eau:kg / (m**2 s)
77 REAL dflux_s(klon) ! derivee du flux sensible dF / dTs
78 REAL dflux_l(klon) ! derivee du flux latent dF / dTs
79
80 ! Flux d'eau "perdue" par la surface et n\'ecessaire pour que limiter la
81 ! hauteur de neige, en kg / m2 / s
82 REAL fqcalving(klon)
83
84 ! Flux thermique utiliser pour fondre la neige
85 REAL ffonte(klon)
86
87 REAL run_off_lic_0(klon)! runof glacier au pas de temps precedent
88
89 ! Local:
90
91 REAL evap(klon) ! evaporation au sol
92
93 INTEGER i, k
94 REAL zx_cq(klon, klev)
95 REAL zx_dq(klon, klev)
96 REAL zx_ch(klon, klev)
97 REAL zx_dh(klon, klev)
98 REAL zx_buf1(klon)
99 REAL zx_buf2(klon)
100 REAL zx_coef(klon, klev)
101 REAL local_h(klon, klev) ! enthalpie potentielle
102 REAL local_q(klon, klev)
103 REAL psref(klon) ! pression de reference pour temperature potent.
104 REAL zx_pkh(klon, klev), zx_pkf(klon, klev)
105
106 ! contre-gradient pour la vapeur d'eau: (kg / kg) / metre
107 REAL gamq(klon, 2:klev)
108 ! contre-gradient pour la chaleur sensible: Kelvin / metre
109 REAL gamt(klon, 2:klev)
110 REAL z_gamaq(klon, 2:klev), z_gamah(klon, 2:klev)
111 REAL zdelz
112
113 real temp_air(klon), spechum(klon)
114 real tq_cdrag(klon), petAcoef(klon), peqAcoef(klon)
115 real petBcoef(klon), peqBcoef(klon)
116 real p1lay(klon)
117
118 real fluxsens(klon)
119 real tsurf_new(knon)
120 real zzpk
121
122 !----------------------------------------------------------------
123
124 if (iflag_pbl == 1) then
125 do k = 3, klev
126 do i = 1, knon
127 gamq(i, k)= 0.0
128 gamt(i, k)= - 1.0e-03
129 enddo
130 enddo
131 do i = 1, knon
132 gamq(i, 2) = 0.0
133 gamt(i, 2) = - 2.5e-03
134 enddo
135 else
136 do k = 2, klev
137 do i = 1, knon
138 gamq(i, k) = 0.0
139 gamt(i, k) = 0.0
140 enddo
141 enddo
142 endif
143
144 DO i = 1, knon
145 psref(i) = paprs(i, 1) !pression de reference est celle au sol
146 ENDDO
147 DO k = 1, klev
148 DO i = 1, knon
149 zx_pkh(i, k) = (psref(i) / paprs(i, k))**RKAPPA
150 zx_pkf(i, k) = (psref(i) / pplay(i, k))**RKAPPA
151 local_h(i, k) = RCPD * t(i, k) * zx_pkf(i, k)
152 local_q(i, k) = q(i, k)
153 ENDDO
154 ENDDO
155
156 ! Convertir les coefficients en variables convenables au calcul:
157
158 DO k = 2, klev
159 DO i = 1, knon
160 zx_coef(i, k) = coef(i, k)*RG / (pplay(i, k - 1) - pplay(i, k)) &
161 *(paprs(i, k)*2 / (t(i, k)+t(i, k - 1)) / RD)**2
162 zx_coef(i, k) = zx_coef(i, k) * dtime*RG
163 ENDDO
164 ENDDO
165
166 ! Preparer les flux lies aux contre-gardients
167
168 DO k = 2, klev
169 DO i = 1, knon
170 zdelz = RD * (t(i, k - 1)+t(i, k)) / 2.0 / RG / paprs(i, k) &
171 *(pplay(i, k - 1) - pplay(i, k))
172 z_gamaq(i, k) = gamq(i, k) * zdelz
173 z_gamah(i, k) = gamt(i, k) * zdelz *RCPD * zx_pkh(i, k)
174 ENDDO
175 ENDDO
176 DO i = 1, knon
177 zx_buf1(i) = zx_coef(i, klev) + delp(i, klev)
178 zx_cq(i, klev) = (local_q(i, klev)*delp(i, klev) &
179 - zx_coef(i, klev)*z_gamaq(i, klev)) / zx_buf1(i)
180 zx_dq(i, klev) = zx_coef(i, klev) / zx_buf1(i)
181
182 zzpk=(pplay(i, klev) / psref(i))**RKAPPA
183 zx_buf2(i) = zzpk*delp(i, klev) + zx_coef(i, klev)
184 zx_ch(i, klev) = (local_h(i, klev)*zzpk*delp(i, klev) &
185 - zx_coef(i, klev)*z_gamah(i, klev)) / zx_buf2(i)
186 zx_dh(i, klev) = zx_coef(i, klev) / zx_buf2(i)
187 ENDDO
188 DO k = klev - 1, 2, - 1
189 DO i = 1, knon
190 zx_buf1(i) = delp(i, k)+zx_coef(i, k) &
191 +zx_coef(i, k+1)*(1. - zx_dq(i, k+1))
192 zx_cq(i, k) = (local_q(i, k)*delp(i, k) &
193 +zx_coef(i, k+1)*zx_cq(i, k+1) &
194 +zx_coef(i, k+1)*z_gamaq(i, k+1) &
195 - zx_coef(i, k)*z_gamaq(i, k)) / zx_buf1(i)
196 zx_dq(i, k) = zx_coef(i, k) / zx_buf1(i)
197
198 zzpk=(pplay(i, k) / psref(i))**RKAPPA
199 zx_buf2(i) = zzpk*delp(i, k)+zx_coef(i, k) &
200 +zx_coef(i, k+1)*(1. - zx_dh(i, k+1))
201 zx_ch(i, k) = (local_h(i, k)*zzpk*delp(i, k) &
202 +zx_coef(i, k+1)*zx_ch(i, k+1) &
203 +zx_coef(i, k+1)*z_gamah(i, k+1) &
204 - zx_coef(i, k)*z_gamah(i, k)) / zx_buf2(i)
205 zx_dh(i, k) = zx_coef(i, k) / zx_buf2(i)
206 ENDDO
207 ENDDO
208
209 DO i = 1, knon
210 zx_buf1(i) = delp(i, 1) + zx_coef(i, 2)*(1. - zx_dq(i, 2))
211 zx_cq(i, 1) = (local_q(i, 1)*delp(i, 1) &
212 +zx_coef(i, 2)*(z_gamaq(i, 2)+zx_cq(i, 2))) &
213 / zx_buf1(i)
214 zx_dq(i, 1) = - 1. * RG / zx_buf1(i)
215
216 zzpk=(pplay(i, 1) / psref(i))**RKAPPA
217 zx_buf2(i) = zzpk*delp(i, 1) + zx_coef(i, 2)*(1. - zx_dh(i, 2))
218 zx_ch(i, 1) = (local_h(i, 1)*zzpk*delp(i, 1) &
219 +zx_coef(i, 2)*(z_gamah(i, 2)+zx_ch(i, 2))) &
220 / zx_buf2(i)
221 zx_dh(i, 1) = - 1. * RG / zx_buf2(i)
222 ENDDO
223
224 ! Appel a interfsurf (appel generique) routine d'interface avec la surface
225
226 ! initialisation
227 petAcoef =0.
228 peqAcoef = 0.
229 petBcoef =0.
230 peqBcoef = 0.
231 p1lay =0.
232
233 petAcoef(1:knon) = zx_ch(1:knon, 1)
234 peqAcoef(1:knon) = zx_cq(1:knon, 1)
235 petBcoef(1:knon) = zx_dh(1:knon, 1)
236 peqBcoef(1:knon) = zx_dq(1:knon, 1)
237 tq_cdrag(1:knon) =coef(:knon, 1)
238 temp_air(1:knon) =t(1:knon, 1)
239 spechum(1:knon)=q(1:knon, 1)
240 p1lay(1:knon) = pplay(1:knon, 1)
241
242 CALL interfsurf_hq(dtime, jour, rmu0, nisurf, knon, knindex, rlat, debut, &
243 nsoilmx, tsoil, qsol, u1lay, v1lay, temp_air, spechum, tq_cdrag, &
244 petAcoef, peqAcoef, petBcoef, peqBcoef, precip_rain, precip_snow, &
245 fder, rugos, rugoro, snow, qsurf, ts(:knon), p1lay, psref, radsol, &
246 evap, fluxsens, fluxlat, dflux_l, dflux_s, tsurf_new, albedo, &
247 z0_new, pctsrf_new_sic, agesno, fqcalving, ffonte, run_off_lic_0)
248
249 flux_t(:knon, 1) = fluxsens(:knon)
250 flux_q(:knon, 1) = - evap(:knon)
251 d_ts = tsurf_new - ts(:knon)
252
253 !==== une fois on a zx_h_ts, on peut faire l'iteration ========
254 DO i = 1, knon
255 local_h(i, 1) = zx_ch(i, 1) + zx_dh(i, 1)*flux_t(i, 1)*dtime
256 local_q(i, 1) = zx_cq(i, 1) + zx_dq(i, 1)*flux_q(i, 1)*dtime
257 ENDDO
258 DO k = 2, klev
259 DO i = 1, knon
260 local_q(i, k) = zx_cq(i, k) + zx_dq(i, k)*local_q(i, k - 1)
261 local_h(i, k) = zx_ch(i, k) + zx_dh(i, k)*local_h(i, k - 1)
262 ENDDO
263 ENDDO
264
265 !== flux_q est le flux de vapeur d'eau: kg / (m**2 s) positive vers bas
266 !== flux_t est le flux de cpt (energie sensible): j / (m**2 s)
267 DO k = 2, klev
268 DO i = 1, knon
269 flux_q(i, k) = (zx_coef(i, k) / RG / dtime) &
270 * (local_q(i, k) - local_q(i, k - 1)+z_gamaq(i, k))
271 flux_t(i, k) = (zx_coef(i, k) / RG / dtime) &
272 * (local_h(i, k) - local_h(i, k - 1)+z_gamah(i, k)) &
273 / zx_pkh(i, k)
274 ENDDO
275 ENDDO
276
277 ! Calcul tendances
278 DO k = 1, klev
279 DO i = 1, knon
280 d_t(i, k) = local_h(i, k) / zx_pkf(i, k) / RCPD - t(i, k)
281 d_q(i, k) = local_q(i, k) - q(i, k)
282 ENDDO
283 ENDDO
284
285 END SUBROUTINE clqh
286
287 end module clqh_m
  ViewVC Help
Powered by ViewVC 1.1.21