/[lmdze]/trunk/Sources/phylmd/clmain.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/clmain.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 40 - (show annotations)
Tue Feb 22 13:49:36 2011 UTC (13 years, 1 month ago) by guez
Original Path: trunk/libf/phylmd/clmain.f90
File size: 27745 byte(s)
"alpha" useless, always 0, in "exner_hyb".

1 module clmain_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE clmain(dtime, itap, date0, pctsrf, pctsrf_new, t, q, u, v,&
8 jour, rmu0, co2_ppm, ok_veget, ocean, npas, nexca, ts,&
9 soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil,&
10 qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat,&
11 rain_f, snow_f, solsw, sollw, sollwdown, fder, rlon, rlat, cufi,&
12 cvfi, rugos, debut, lafin, agesno, rugoro, d_t, d_q, d_u, d_v,&
13 d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2,&
14 dflux_t, dflux_q, zcoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh,&
15 capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl,&
16 fqcalving, ffonte, run_off_lic_0, flux_o, flux_g, tslab, seaice)
17
18 ! From phylmd/clmain.F, version 1.6 2005/11/16 14:47:19
19 ! Author: Z.X. Li (LMD/CNRS), date: 1993/08/18
20 ! Objet : interface de "couche limite" (diffusion verticale)
21
22 ! Tout ce qui a trait aux traceurs est dans "phytrac" maintenant.
23 ! Pour l'instant le calcul de la couche limite pour les traceurs
24 ! se fait avec "cltrac" et ne tient pas compte de la différentiation
25 ! des sous-fractions de sol.
26
27 ! Pour pouvoir extraire les coefficients d'échanges et le vent
28 ! dans la première couche, trois champs supplémentaires ont été
29 ! créés : "zcoefh", "zu1" et "zv1". Pour l'instant nous avons
30 ! moyenné les valeurs de ces trois champs sur les 4 sous-surfaces
31 ! du modèle. Dans l'avenir, si les informations des sous-surfaces
32 ! doivent être prises en compte, il faudra sortir ces mêmes champs
33 ! en leur ajoutant une dimension, c'est-à-dire "nbsrf" (nombre de
34 ! sous-surfaces).
35
36 ! Arguments:
37 ! dtime----input-R- interval du temps (secondes)
38 ! itap-----input-I- numero du pas de temps
39 ! date0----input-R- jour initial
40 ! t--------input-R- temperature (K)
41 ! q--------input-R- vapeur d'eau (kg/kg)
42 ! u--------input-R- vitesse u
43 ! v--------input-R- vitesse v
44 ! ts-------input-R- temperature du sol (en Kelvin)
45 ! paprs----input-R- pression a intercouche (Pa)
46 ! pplay----input-R- pression au milieu de couche (Pa)
47 ! radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2
48 ! rlat-----input-R- latitude en degree
49 ! rugos----input-R- longeur de rugosite (en m)
50 ! cufi-----input-R- resolution des mailles en x (m)
51 ! cvfi-----input-R- resolution des mailles en y (m)
52
53 ! d_t------output-R- le changement pour "t"
54 ! d_q------output-R- le changement pour "q"
55 ! d_u------output-R- le changement pour "u"
56 ! d_v------output-R- le changement pour "v"
57 ! d_ts-----output-R- le changement pour "ts"
58 ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
59 ! (orientation positive vers le bas)
60 ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
61 ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
62 ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
63 ! dflux_t derive du flux sensible
64 ! dflux_q derive du flux latent
65 !IM "slab" ocean
66 ! flux_g---output-R- flux glace (pour OCEAN='slab ')
67 ! flux_o---output-R- flux ocean (pour OCEAN='slab ')
68
69 ! tslab-in/output-R temperature du slab ocean (en Kelvin)
70 ! uniqmnt pour slab
71
72 ! seaice---output-R- glace de mer (kg/m2) (pour OCEAN='slab ')
73 !cc
74 ! ffonte----Flux thermique utilise pour fondre la neige
75 ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
76 ! hauteur de neige, en kg/m2/s
77 ! on rajoute en output yu1 et yv1 qui sont les vents dans
78 ! la premiere couche
79 ! ces 4 variables sont maintenant traites dans phytrac
80 ! itr--------input-I- nombre de traceurs
81 ! tr---------input-R- q. de traceurs
82 ! flux_surf--input-R- flux de traceurs a la surface
83 ! d_tr-------output-R tendance de traceurs
84 !IM cf. AM : PBL
85 ! trmb1-------deep_cape
86 ! trmb2--------inhibition
87 ! trmb3-------Point Omega
88 ! Cape(klon)-------Cape du thermique
89 ! EauLiq(klon)-------Eau liqu integr du thermique
90 ! ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL
91 ! lcl------- Niveau de condensation
92 ! pblh------- HCL
93 ! pblT------- T au nveau HCL
94
95 USE histcom, ONLY : histbeg_totreg, histdef, histend, histsync
96 use histwrite_m, only: histwrite
97 use calendar, ONLY : ymds2ju
98 USE dimens_m, ONLY : iim, jjm
99 USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
100 USE dimphy, ONLY : klev, klon, zmasq
101 USE dimsoil, ONLY : nsoilmx
102 USE temps, ONLY : annee_ref, itau_phy
103 USE dynetat0_m, ONLY : day_ini
104 USE iniprint, ONLY : prt_level
105 USE suphec_m, ONLY : rd, rg, rkappa
106 USE conf_phys_m, ONLY : iflag_pbl
107 USE gath_cpl, ONLY : gath2cpl
108 use hbtm_m, only: hbtm
109
110 REAL, INTENT (IN) :: dtime
111 REAL date0
112 INTEGER, INTENT (IN) :: itap
113 REAL t(klon, klev), q(klon, klev)
114 REAL u(klon, klev), v(klon, klev)
115 REAL, INTENT (IN) :: paprs(klon, klev+1)
116 REAL, INTENT (IN) :: pplay(klon, klev)
117 REAL, INTENT (IN) :: rlon(klon), rlat(klon)
118 REAL cufi(klon), cvfi(klon)
119 REAL d_t(klon, klev), d_q(klon, klev)
120 REAL d_u(klon, klev), d_v(klon, klev)
121 REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
122 REAL dflux_t(klon), dflux_q(klon)
123 !IM "slab" ocean
124 REAL flux_o(klon), flux_g(klon)
125 REAL y_flux_o(klon), y_flux_g(klon)
126 REAL tslab(klon), ytslab(klon)
127 REAL seaice(klon), y_seaice(klon)
128 REAL y_fqcalving(klon), y_ffonte(klon)
129 REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
130 REAL run_off_lic_0(klon), y_run_off_lic_0(klon)
131
132 REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
133 REAL rugmer(klon), agesno(klon, nbsrf)
134 REAL, INTENT (IN) :: rugoro(klon)
135 REAL cdragh(klon), cdragm(klon)
136 ! jour de l'annee en cours
137 INTEGER jour
138 REAL rmu0(klon) ! cosinus de l'angle solaire zenithal
139 ! taux CO2 atmosphere
140 REAL co2_ppm
141 LOGICAL, INTENT (IN) :: debut
142 LOGICAL, INTENT (IN) :: lafin
143 LOGICAL ok_veget
144 CHARACTER (len=*), INTENT (IN) :: ocean
145 INTEGER npas, nexca
146
147 REAL pctsrf(klon, nbsrf)
148 REAL ts(klon, nbsrf)
149 REAL d_ts(klon, nbsrf)
150 REAL snow(klon, nbsrf)
151 REAL qsurf(klon, nbsrf)
152 REAL evap(klon, nbsrf)
153 REAL albe(klon, nbsrf)
154 REAL alblw(klon, nbsrf)
155
156 REAL fluxlat(klon, nbsrf)
157
158 REAL rain_f(klon), snow_f(klon)
159 REAL fder(klon)
160
161 REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)
162 REAL rugos(klon, nbsrf)
163 ! la nouvelle repartition des surfaces sortie de l'interface
164 REAL pctsrf_new(klon, nbsrf)
165
166 REAL zcoefh(klon, klev)
167 REAL zu1(klon)
168 REAL zv1(klon)
169
170 !$$$ PB ajout pour soil
171 LOGICAL, INTENT (IN) :: soil_model
172 !IM ajout seuils cdrm, cdrh
173 REAL cdmmax, cdhmax
174
175 REAL ksta, ksta_ter
176 LOGICAL ok_kzmin
177
178 REAL ftsoil(klon, nsoilmx, nbsrf)
179 REAL ytsoil(klon, nsoilmx)
180 REAL qsol(klon)
181
182 EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
183
184 REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
185 REAL yalb(klon)
186 REAL yalblw(klon)
187 REAL yu1(klon), yv1(klon)
188 REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
189 REAL yrain_f(klon), ysnow_f(klon)
190 REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)
191 REAL yfder(klon), ytaux(klon), ytauy(klon)
192 REAL yrugm(klon), yrads(klon), yrugoro(klon)
193
194 REAL yfluxlat(klon)
195
196 REAL y_d_ts(klon)
197 REAL y_d_t(klon, klev), y_d_q(klon, klev)
198 REAL y_d_u(klon, klev), y_d_v(klon, klev)
199 REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
200 REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
201 REAL y_dflux_t(klon), y_dflux_q(klon)
202 REAL ycoefh(klon, klev), ycoefm(klon, klev)
203 REAL yu(klon, klev), yv(klon, klev)
204 REAL yt(klon, klev), yq(klon, klev)
205 REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
206
207 LOGICAL ok_nonloc
208 PARAMETER (ok_nonloc=.FALSE.)
209 REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
210
211 !IM 081204 hcl_Anne ? BEG
212 REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
213 REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
214 REAL ykmq(klon, klev+1)
215 REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)
216 REAL q2diag(klon, klev+1)
217 !IM 081204 hcl_Anne ? END
218
219 REAL u1lay(klon), v1lay(klon)
220 REAL delp(klon, klev)
221 INTEGER i, k, nsrf
222
223 INTEGER ni(klon), knon, j
224
225 REAL pctsrf_pot(klon, nbsrf)
226 ! "pourcentage potentiel" pour tenir compte des éventuelles
227 ! apparitions ou disparitions de la glace de mer
228
229 REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
230
231 ! maf pour sorties IOISPL en cas de debugagage
232
233 CHARACTER (80) cldebug
234 SAVE cldebug
235 CHARACTER (8) cl_surf(nbsrf)
236 SAVE cl_surf
237 INTEGER nhoridbg, nidbg
238 SAVE nhoridbg, nidbg
239 INTEGER ndexbg(iim*(jjm+1))
240 REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian
241 REAL tabindx(klon)
242 REAL debugtab(iim, jjm+1)
243 LOGICAL first_appel
244 SAVE first_appel
245 DATA first_appel/ .TRUE./
246 LOGICAL :: debugindex = .FALSE.
247 INTEGER idayref
248 REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
249 REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
250
251 REAL yt2m(klon), yq2m(klon), yu10m(klon)
252 REAL yustar(klon)
253 ! -- LOOP
254 REAL yu10mx(klon)
255 REAL yu10my(klon)
256 REAL ywindsp(klon)
257 ! -- LOOP
258
259 REAL yt10m(klon), yq10m(klon)
260 !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
261 ! physiq ce qui permet de sortir les grdeurs par sous surface)
262 REAL pblh(klon, nbsrf)
263 REAL plcl(klon, nbsrf)
264 REAL capcl(klon, nbsrf)
265 REAL oliqcl(klon, nbsrf)
266 REAL cteicl(klon, nbsrf)
267 REAL pblt(klon, nbsrf)
268 REAL therm(klon, nbsrf)
269 REAL trmb1(klon, nbsrf)
270 REAL trmb2(klon, nbsrf)
271 REAL trmb3(klon, nbsrf)
272 REAL ypblh(klon)
273 REAL ylcl(klon)
274 REAL ycapcl(klon)
275 REAL yoliqcl(klon)
276 REAL ycteicl(klon)
277 REAL ypblt(klon)
278 REAL ytherm(klon)
279 REAL ytrmb1(klon)
280 REAL ytrmb2(klon)
281 REAL ytrmb3(klon)
282 REAL y_cd_h(klon), y_cd_m(klon)
283 REAL uzon(klon), vmer(klon)
284 REAL tair1(klon), qair1(klon), tairsol(klon)
285 REAL psfce(klon), patm(klon)
286
287 REAL qairsol(klon), zgeo1(klon)
288 REAL rugo1(klon)
289
290 ! utiliser un jeu de fonctions simples
291 LOGICAL zxli
292 PARAMETER (zxli=.FALSE.)
293
294 REAL zt, zqs, zdelta, zcor
295 REAL t_coup
296 PARAMETER (t_coup=273.15)
297
298 CHARACTER (len=20) :: modname = 'clmain'
299
300 !------------------------------------------------------------
301
302 ytherm = 0.
303
304 IF (debugindex .AND. first_appel) THEN
305 first_appel = .FALSE.
306
307 ! initialisation sorties netcdf
308
309 idayref = day_ini
310 CALL ymds2ju(annee_ref, 1, idayref, 0., zjulian)
311 CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)
312 DO i = 1, iim
313 zx_lon(i, 1) = rlon(i+1)
314 zx_lon(i, jjm+1) = rlon(i+1)
315 END DO
316 CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlat, zx_lat)
317 cldebug = 'sous_index'
318 CALL histbeg_totreg(cldebug, zx_lon(:, 1), zx_lat(1, :), 1, &
319 iim, 1, jjm+1, itau_phy, zjulian, dtime, nhoridbg, nidbg)
320 ! no vertical axis
321 cl_surf(1) = 'ter'
322 cl_surf(2) = 'lic'
323 cl_surf(3) = 'oce'
324 cl_surf(4) = 'sic'
325 DO nsrf = 1, nbsrf
326 CALL histdef(nidbg, cl_surf(nsrf), cl_surf(nsrf), '-', iim, jjm+1, &
327 nhoridbg, 1, 1, 1, -99, 'inst', dtime, dtime)
328 END DO
329 CALL histend(nidbg)
330 CALL histsync(nidbg)
331 END IF
332
333 DO k = 1, klev ! epaisseur de couche
334 DO i = 1, klon
335 delp(i, k) = paprs(i, k) - paprs(i, k+1)
336 END DO
337 END DO
338 DO i = 1, klon ! vent de la premiere couche
339 zx_alf1 = 1.0
340 zx_alf2 = 1.0 - zx_alf1
341 u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2
342 v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
343 END DO
344
345 ! Initialization:
346 rugmer = 0.
347 cdragh = 0.
348 cdragm = 0.
349 dflux_t = 0.
350 dflux_q = 0.
351 zu1 = 0.
352 zv1 = 0.
353 ypct = 0.
354 yts = 0.
355 ysnow = 0.
356 yqsurf = 0.
357 yalb = 0.
358 yalblw = 0.
359 yrain_f = 0.
360 ysnow_f = 0.
361 yfder = 0.
362 ytaux = 0.
363 ytauy = 0.
364 ysolsw = 0.
365 ysollw = 0.
366 ysollwdown = 0.
367 yrugos = 0.
368 yu1 = 0.
369 yv1 = 0.
370 yrads = 0.
371 ypaprs = 0.
372 ypplay = 0.
373 ydelp = 0.
374 yu = 0.
375 yv = 0.
376 yt = 0.
377 yq = 0.
378 pctsrf_new = 0.
379 y_flux_u = 0.
380 y_flux_v = 0.
381 !$$ PB
382 y_dflux_t = 0.
383 y_dflux_q = 0.
384 ytsoil = 999999.
385 yrugoro = 0.
386 ! -- LOOP
387 yu10mx = 0.
388 yu10my = 0.
389 ywindsp = 0.
390 ! -- LOOP
391 d_ts = 0.
392 !§§§ PB
393 yfluxlat = 0.
394 flux_t = 0.
395 flux_q = 0.
396 flux_u = 0.
397 flux_v = 0.
398 d_t = 0.
399 d_q = 0.
400 d_u = 0.
401 d_v = 0.
402 zcoefh = 0.
403
404 ! Boucler sur toutes les sous-fractions du sol:
405
406 ! Initialisation des "pourcentages potentiels". On considère ici qu'on
407 ! peut avoir potentiellement de la glace sur tout le domaine océanique
408 ! (à affiner)
409
410 pctsrf_pot = pctsrf
411 pctsrf_pot(:, is_oce) = 1. - zmasq
412 pctsrf_pot(:, is_sic) = 1. - zmasq
413
414 DO nsrf = 1, nbsrf
415 ! chercher les indices:
416 ni = 0
417 knon = 0
418 DO i = 1, klon
419 ! Pour déterminer le domaine à traiter, on utilise les surfaces
420 ! "potentielles"
421 IF (pctsrf_pot(i, nsrf) > epsfra) THEN
422 knon = knon + 1
423 ni(knon) = i
424 END IF
425 END DO
426
427 ! variables pour avoir une sortie IOIPSL des INDEX
428 IF (debugindex) THEN
429 tabindx = 0.
430 DO i = 1, knon
431 tabindx(i) = real(i)
432 END DO
433 debugtab = 0.
434 ndexbg = 0
435 CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)
436 CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
437 END IF
438
439 IF (knon==0) CYCLE
440
441 DO j = 1, knon
442 i = ni(j)
443 ypct(j) = pctsrf(i, nsrf)
444 yts(j) = ts(i, nsrf)
445 ytslab(i) = tslab(i)
446 ysnow(j) = snow(i, nsrf)
447 yqsurf(j) = qsurf(i, nsrf)
448 yalb(j) = albe(i, nsrf)
449 yalblw(j) = alblw(i, nsrf)
450 yrain_f(j) = rain_f(i)
451 ysnow_f(j) = snow_f(i)
452 yagesno(j) = agesno(i, nsrf)
453 yfder(j) = fder(i)
454 ytaux(j) = flux_u(i, 1, nsrf)
455 ytauy(j) = flux_v(i, 1, nsrf)
456 ysolsw(j) = solsw(i, nsrf)
457 ysollw(j) = sollw(i, nsrf)
458 ysollwdown(j) = sollwdown(i)
459 yrugos(j) = rugos(i, nsrf)
460 yrugoro(j) = rugoro(i)
461 yu1(j) = u1lay(i)
462 yv1(j) = v1lay(i)
463 yrads(j) = ysolsw(j) + ysollw(j)
464 ypaprs(j, klev+1) = paprs(i, klev+1)
465 y_run_off_lic_0(j) = run_off_lic_0(i)
466 yu10mx(j) = u10m(i, nsrf)
467 yu10my(j) = v10m(i, nsrf)
468 ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
469 END DO
470
471 ! IF bucket model for continent, copy soil water content
472 IF (nsrf==is_ter .AND. .NOT. ok_veget) THEN
473 DO j = 1, knon
474 i = ni(j)
475 yqsol(j) = qsol(i)
476 END DO
477 ELSE
478 yqsol = 0.
479 END IF
480 !$$$ PB ajour pour soil
481 DO k = 1, nsoilmx
482 DO j = 1, knon
483 i = ni(j)
484 ytsoil(j, k) = ftsoil(i, k, nsrf)
485 END DO
486 END DO
487 DO k = 1, klev
488 DO j = 1, knon
489 i = ni(j)
490 ypaprs(j, k) = paprs(i, k)
491 ypplay(j, k) = pplay(i, k)
492 ydelp(j, k) = delp(i, k)
493 yu(j, k) = u(i, k)
494 yv(j, k) = v(i, k)
495 yt(j, k) = t(i, k)
496 yq(j, k) = q(i, k)
497 END DO
498 END DO
499
500 ! calculer Cdrag et les coefficients d'echange
501 CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&
502 yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)
503 !IM 081204 BEG
504 !CR test
505 IF (iflag_pbl==1) THEN
506 !IM 081204 END
507 CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
508 DO k = 1, klev
509 DO i = 1, knon
510 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
511 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
512 END DO
513 END DO
514 END IF
515
516 !IM cf JLD : on seuille ycoefm et ycoefh
517 IF (nsrf==is_oce) THEN
518 DO j = 1, knon
519 ! ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)
520 ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)
521 ! ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)
522 ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)
523 END DO
524 END IF
525
526 !IM: 261103
527 IF (ok_kzmin) THEN
528 !IM cf FH: 201103 BEG
529 ! Calcul d'une diffusion minimale pour les conditions tres stables.
530 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, &
531 ycoefm0, ycoefh0)
532
533 IF (1==1) THEN
534 DO k = 1, klev
535 DO i = 1, knon
536 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
537 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
538 END DO
539 END DO
540 END IF
541 !IM cf FH: 201103 END
542 !IM: 261103
543 END IF !ok_kzmin
544
545 IF (iflag_pbl>=3) THEN
546 ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin
547 yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &
548 1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg
549 DO k = 2, klev
550 yzlay(1:knon, k) = yzlay(1:knon, k-1) &
551 + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
552 / ypaprs(1:knon, k) &
553 * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
554 END DO
555 DO k = 1, klev
556 yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
557 / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
558 END DO
559 yzlev(1:knon, 1) = 0.
560 yzlev(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1:knon, klev-1)
561 DO k = 2, klev
562 yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
563 END DO
564 DO k = 1, klev + 1
565 DO j = 1, knon
566 i = ni(j)
567 yq2(j, k) = q2(i, k, nsrf)
568 END DO
569 END DO
570
571 ! Bug introduit volontairement pour converger avec les resultats
572 ! du papier sur les thermiques.
573 IF (1==1) THEN
574 y_cd_m(1:knon) = ycoefm(1:knon, 1)
575 y_cd_h(1:knon) = ycoefh(1:knon, 1)
576 ELSE
577 y_cd_h(1:knon) = ycoefm(1:knon, 1)
578 y_cd_m(1:knon) = ycoefh(1:knon, 1)
579 END IF
580 CALL ustarhb(knon, yu, yv, y_cd_m, yustar)
581
582 IF (prt_level>9) THEN
583 PRINT *, 'USTAR = ', yustar
584 END IF
585
586 ! iflag_pbl peut etre utilise comme longuer de melange
587
588 IF (iflag_pbl>=11) THEN
589 CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
590 yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &
591 iflag_pbl)
592 ELSE
593 CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &
594 yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
595 END IF
596
597 ycoefm(1:knon, 1) = y_cd_m(1:knon)
598 ycoefh(1:knon, 1) = y_cd_h(1:knon)
599 ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)
600 ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)
601 END IF
602
603 ! calculer la diffusion des vitesses "u" et "v"
604 CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &
605 ydelp, y_d_u, y_flux_u)
606 CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &
607 ydelp, y_d_v, y_flux_v)
608
609 ! pour le couplage
610 ytaux = y_flux_u(:, 1)
611 ytauy = y_flux_v(:, 1)
612
613 ! calculer la diffusion de "q" et de "h"
614 CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&
615 cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&
616 yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&
617 yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&
618 ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &
619 yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&
620 yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&
621 yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&
622 y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&
623 ytslab, y_seaice)
624
625 ! calculer la longueur de rugosite sur ocean
626 yrugm = 0.
627 IF (nsrf==is_oce) THEN
628 DO j = 1, knon
629 yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
630 0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))
631 yrugm(j) = max(1.5E-05, yrugm(j))
632 END DO
633 END IF
634 DO j = 1, knon
635 y_dflux_t(j) = y_dflux_t(j)*ypct(j)
636 y_dflux_q(j) = y_dflux_q(j)*ypct(j)
637 yu1(j) = yu1(j)*ypct(j)
638 yv1(j) = yv1(j)*ypct(j)
639 END DO
640
641 DO k = 1, klev
642 DO j = 1, knon
643 i = ni(j)
644 ycoefh(j, k) = ycoefh(j, k)*ypct(j)
645 ycoefm(j, k) = ycoefm(j, k)*ypct(j)
646 y_d_t(j, k) = y_d_t(j, k)*ypct(j)
647 y_d_q(j, k) = y_d_q(j, k)*ypct(j)
648 !§§§ PB
649 flux_t(i, k, nsrf) = y_flux_t(j, k)
650 flux_q(i, k, nsrf) = y_flux_q(j, k)
651 flux_u(i, k, nsrf) = y_flux_u(j, k)
652 flux_v(i, k, nsrf) = y_flux_v(j, k)
653 !$$$ PB y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)
654 !$$$ PB y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)
655 y_d_u(j, k) = y_d_u(j, k)*ypct(j)
656 y_d_v(j, k) = y_d_v(j, k)*ypct(j)
657 !$$$ PB y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)
658 !$$$ PB y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)
659 END DO
660 END DO
661
662 evap(:, nsrf) = -flux_q(:, 1, nsrf)
663
664 albe(:, nsrf) = 0.
665 alblw(:, nsrf) = 0.
666 snow(:, nsrf) = 0.
667 qsurf(:, nsrf) = 0.
668 rugos(:, nsrf) = 0.
669 fluxlat(:, nsrf) = 0.
670 DO j = 1, knon
671 i = ni(j)
672 d_ts(i, nsrf) = y_d_ts(j)
673 albe(i, nsrf) = yalb(j)
674 alblw(i, nsrf) = yalblw(j)
675 snow(i, nsrf) = ysnow(j)
676 qsurf(i, nsrf) = yqsurf(j)
677 rugos(i, nsrf) = yz0_new(j)
678 fluxlat(i, nsrf) = yfluxlat(j)
679 !$$$ pb rugmer(i) = yrugm(j)
680 IF (nsrf==is_oce) THEN
681 rugmer(i) = yrugm(j)
682 rugos(i, nsrf) = yrugm(j)
683 END IF
684 !IM cf JLD ??
685 agesno(i, nsrf) = yagesno(j)
686 fqcalving(i, nsrf) = y_fqcalving(j)
687 ffonte(i, nsrf) = y_ffonte(j)
688 cdragh(i) = cdragh(i) + ycoefh(j, 1)
689 cdragm(i) = cdragm(i) + ycoefm(j, 1)
690 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
691 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
692 zu1(i) = zu1(i) + yu1(j)
693 zv1(i) = zv1(i) + yv1(j)
694 END DO
695 IF (nsrf==is_ter) THEN
696 DO j = 1, knon
697 i = ni(j)
698 qsol(i) = yqsol(j)
699 END DO
700 END IF
701 IF (nsrf==is_lic) THEN
702 DO j = 1, knon
703 i = ni(j)
704 run_off_lic_0(i) = y_run_off_lic_0(j)
705 END DO
706 END IF
707 !$$$ PB ajout pour soil
708 ftsoil(:, :, nsrf) = 0.
709 DO k = 1, nsoilmx
710 DO j = 1, knon
711 i = ni(j)
712 ftsoil(i, k, nsrf) = ytsoil(j, k)
713 END DO
714 END DO
715
716 DO j = 1, knon
717 i = ni(j)
718 DO k = 1, klev
719 d_t(i, k) = d_t(i, k) + y_d_t(j, k)
720 d_q(i, k) = d_q(i, k) + y_d_q(j, k)
721 !$$$ PB flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)
722 !$$$ flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k)
723 d_u(i, k) = d_u(i, k) + y_d_u(j, k)
724 d_v(i, k) = d_v(i, k) + y_d_v(j, k)
725 !$$$ PB flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)
726 !$$$ flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)
727 zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)
728 END DO
729 END DO
730
731 !cc diagnostic t, q a 2m et u, v a 10m
732
733 DO j = 1, knon
734 i = ni(j)
735 uzon(j) = yu(j, 1) + y_d_u(j, 1)
736 vmer(j) = yv(j, 1) + y_d_v(j, 1)
737 tair1(j) = yt(j, 1) + y_d_t(j, 1)
738 qair1(j) = yq(j, 1) + y_d_q(j, 1)
739 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
740 1)))*(ypaprs(j, 1)-ypplay(j, 1))
741 tairsol(j) = yts(j) + y_d_ts(j)
742 rugo1(j) = yrugos(j)
743 IF (nsrf==is_oce) THEN
744 rugo1(j) = rugos(i, nsrf)
745 END IF
746 psfce(j) = ypaprs(j, 1)
747 patm(j) = ypplay(j, 1)
748
749 qairsol(j) = yqsurf(j)
750 END DO
751
752 CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &
753 tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &
754 yu10m, yustar)
755 !IM 081204 END
756
757 DO j = 1, knon
758 i = ni(j)
759 t2m(i, nsrf) = yt2m(j)
760 q2m(i, nsrf) = yq2m(j)
761
762 ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
763 u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
764 v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
765
766 END DO
767
768 DO i = 1, knon
769 y_cd_h(i) = ycoefh(i, 1)
770 y_cd_m(i) = ycoefm(i, 1)
771 END DO
772 CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
773 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
774 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
775
776 DO j = 1, knon
777 i = ni(j)
778 pblh(i, nsrf) = ypblh(j)
779 plcl(i, nsrf) = ylcl(j)
780 capcl(i, nsrf) = ycapcl(j)
781 oliqcl(i, nsrf) = yoliqcl(j)
782 cteicl(i, nsrf) = ycteicl(j)
783 pblt(i, nsrf) = ypblt(j)
784 therm(i, nsrf) = ytherm(j)
785 trmb1(i, nsrf) = ytrmb1(j)
786 trmb2(i, nsrf) = ytrmb2(j)
787 trmb3(i, nsrf) = ytrmb3(j)
788 END DO
789
790 DO j = 1, knon
791 DO k = 1, klev + 1
792 i = ni(j)
793 q2(i, k, nsrf) = yq2(j, k)
794 END DO
795 END DO
796 !IM "slab" ocean
797 IF (nsrf==is_oce) THEN
798 DO j = 1, knon
799 ! on projette sur la grille globale
800 i = ni(j)
801 IF (pctsrf_new(i, is_oce)>epsfra) THEN
802 flux_o(i) = y_flux_o(j)
803 ELSE
804 flux_o(i) = 0.
805 END IF
806 END DO
807 END IF
808
809 IF (nsrf==is_sic) THEN
810 DO j = 1, knon
811 i = ni(j)
812 ! On pondère lorsque l'on fait le bilan au sol :
813 ! flux_g(i) = y_flux_g(j)*ypct(j)
814 IF (pctsrf_new(i, is_sic)>epsfra) THEN
815 flux_g(i) = y_flux_g(j)
816 ELSE
817 flux_g(i) = 0.
818 END IF
819 END DO
820
821 END IF
822 !nsrf.EQ.is_sic
823 IF (ocean=='slab ') THEN
824 IF (nsrf==is_oce) THEN
825 tslab(1:klon) = ytslab(1:klon)
826 seaice(1:klon) = y_seaice(1:klon)
827 !nsrf
828 END IF
829 !OCEAN
830 END IF
831 END DO
832
833 ! On utilise les nouvelles surfaces
834 ! A rajouter: conservation de l'albedo
835
836 rugos(:, is_oce) = rugmer
837 pctsrf = pctsrf_new
838
839 END SUBROUTINE clmain
840
841 end module clmain_m

  ViewVC Help
Powered by ViewVC 1.1.21