/[lmdze]/trunk/libf/phylmd/clmain.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/clmain.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 38 - (show annotations)
Thu Jan 6 17:52:19 2011 UTC (13 years, 4 months ago) by guez
File size: 28294 byte(s)
Extracted ASCII art from "inigeom" into a separate text file in the
documentation.

"test_disvert" now creates a separate file for layer thicknesses.

Moved variables from module "yomcst" to module "suphec_m" because this
is where those variables are defined. Kept in "yomcst" only parameters
of Earth orbit. Gave the attribute "parameter" to some variables of
module "suphec_m".

Variables of module "yoethf" were defined in procedure "suphec". Moved
these definitions to a new procedure "yoethf" in module "yoethf_m".

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
20 ! Tout ce qui a trait aux traceurs est dans phytrac maintenant.
21 ! Pour l'instant le calcul de la couche limite pour les traceurs
22 ! se fait avec cltrac et ne tient pas compte de la différentiation
23 ! des sous-fractions de sol.
24
25 ! Pour pouvoir extraire les coefficients d'échanges et le vent
26 ! dans la première couche, trois champs supplémentaires ont été créés :
27 ! zcoefh, zu1 et zv1. Pour l'instant nous avons moyenné les valeurs
28 ! de ces trois champs sur les 4 sous-surfaces du modèle. Dans l'avenir
29 ! si les informations des sous-surfaces doivent être prises en compte
30 ! il faudra sortir ces mêmes champs en leur ajoutant une dimension,
31 ! c'est a dire nbsrf (nombre de sous-surfaces).
32
33 ! Auteur Z.X. Li (LMD/CNRS) date: 1993/08/18
34 ! Objet : interface de "couche limite" (diffusion verticale)
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 ! Introduction d'une variable "pourcentage potentiel" pour tenir compte
225 ! des eventuelles apparitions et/ou disparitions de la glace de mer
226 REAL pctsrf_pot(klon, nbsrf)
227
228 REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
229
230 ! maf pour sorties IOISPL en cas de debugagage
231
232 CHARACTER (80) cldebug
233 SAVE cldebug
234 CHARACTER (8) cl_surf(nbsrf)
235 SAVE cl_surf
236 INTEGER nhoridbg, nidbg
237 SAVE nhoridbg, nidbg
238 INTEGER ndexbg(iim*(jjm+1))
239 REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian
240 REAL tabindx(klon)
241 REAL debugtab(iim, jjm+1)
242 LOGICAL first_appel
243 SAVE first_appel
244 DATA first_appel/ .TRUE./
245 LOGICAL :: debugindex = .FALSE.
246 INTEGER idayref
247 REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
248 REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
249
250 REAL yt2m(klon), yq2m(klon), yu10m(klon)
251 REAL yustar(klon)
252 ! -- LOOP
253 REAL yu10mx(klon)
254 REAL yu10my(klon)
255 REAL ywindsp(klon)
256 ! -- LOOP
257
258 REAL yt10m(klon), yq10m(klon)
259 !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
260 ! physiq ce qui permet de sortir les grdeurs par sous surface)
261 REAL pblh(klon, nbsrf)
262 REAL plcl(klon, nbsrf)
263 REAL capcl(klon, nbsrf)
264 REAL oliqcl(klon, nbsrf)
265 REAL cteicl(klon, nbsrf)
266 REAL pblt(klon, nbsrf)
267 REAL therm(klon, nbsrf)
268 REAL trmb1(klon, nbsrf)
269 REAL trmb2(klon, nbsrf)
270 REAL trmb3(klon, nbsrf)
271 REAL ypblh(klon)
272 REAL ylcl(klon)
273 REAL ycapcl(klon)
274 REAL yoliqcl(klon)
275 REAL ycteicl(klon)
276 REAL ypblt(klon)
277 REAL ytherm(klon)
278 REAL ytrmb1(klon)
279 REAL ytrmb2(klon)
280 REAL ytrmb3(klon)
281 REAL y_cd_h(klon), y_cd_m(klon)
282 REAL uzon(klon), vmer(klon)
283 REAL tair1(klon), qair1(klon), tairsol(klon)
284 REAL psfce(klon), patm(klon)
285
286 REAL qairsol(klon), zgeo1(klon)
287 REAL rugo1(klon)
288
289 ! utiliser un jeu de fonctions simples
290 LOGICAL zxli
291 PARAMETER (zxli=.FALSE.)
292
293 REAL zt, zqs, zdelta, zcor
294 REAL t_coup
295 PARAMETER (t_coup=273.15)
296
297 CHARACTER (len=20) :: modname = 'clmain'
298
299 !------------------------------------------------------------
300
301 ! initialisation Anne
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.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 ! initialisation:
346
347 DO i = 1, klon
348 rugmer(i) = 0.0
349 cdragh(i) = 0.0
350 cdragm(i) = 0.0
351 dflux_t(i) = 0.0
352 dflux_q(i) = 0.0
353 zu1(i) = 0.0
354 zv1(i) = 0.0
355 END DO
356 ypct = 0.0
357 yts = 0.0
358 ysnow = 0.0
359 yqsurf = 0.0
360 yalb = 0.0
361 yalblw = 0.0
362 yrain_f = 0.0
363 ysnow_f = 0.0
364 yfder = 0.0
365 ytaux = 0.0
366 ytauy = 0.0
367 ysolsw = 0.0
368 ysollw = 0.0
369 ysollwdown = 0.0
370 yrugos = 0.0
371 yu1 = 0.0
372 yv1 = 0.0
373 yrads = 0.0
374 ypaprs = 0.0
375 ypplay = 0.0
376 ydelp = 0.0
377 yu = 0.0
378 yv = 0.0
379 yt = 0.0
380 yq = 0.0
381 pctsrf_new = 0.0
382 y_flux_u = 0.0
383 y_flux_v = 0.0
384 !$$ PB
385 y_dflux_t = 0.0
386 y_dflux_q = 0.0
387 ytsoil = 999999.
388 yrugoro = 0.
389 ! -- LOOP
390 yu10mx = 0.0
391 yu10my = 0.0
392 ywindsp = 0.0
393 ! -- LOOP
394 DO nsrf = 1, nbsrf
395 DO i = 1, klon
396 d_ts(i, nsrf) = 0.0
397 END DO
398 END DO
399 !§§§ PB
400 yfluxlat = 0.
401 flux_t = 0.
402 flux_q = 0.
403 flux_u = 0.
404 flux_v = 0.
405 DO k = 1, klev
406 DO i = 1, klon
407 d_t(i, k) = 0.0
408 d_q(i, k) = 0.0
409 d_u(i, k) = 0.0
410 d_v(i, k) = 0.0
411 zcoefh(i, k) = 0.0
412 END DO
413 END DO
414
415 ! Boucler sur toutes les sous-fractions du sol:
416
417 ! Initialisation des "pourcentages potentiels". On considère ici qu'on
418 ! peut avoir potentiellement de la glace sur tout le domaine océanique
419 ! (à affiner)
420
421 pctsrf_pot = pctsrf
422 pctsrf_pot(:, is_oce) = 1. - zmasq
423 pctsrf_pot(:, is_sic) = 1. - zmasq
424
425 DO nsrf = 1, nbsrf
426 ! chercher les indices:
427 ni = 0
428 knon = 0
429 DO i = 1, klon
430 ! pour determiner le domaine a traiter on utilise les surfaces
431 ! "potentielles"
432 IF (pctsrf_pot(i, nsrf) > epsfra) THEN
433 knon = knon + 1
434 ni(knon) = i
435 END IF
436 END DO
437
438 ! variables pour avoir une sortie IOIPSL des INDEX
439 IF (debugindex) THEN
440 tabindx = 0.
441 DO i = 1, knon
442 tabindx(i) = real(i)
443 END DO
444 debugtab = 0.
445 ndexbg = 0
446 CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)
447 CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
448 END IF
449
450 IF (knon==0) CYCLE
451
452 DO j = 1, knon
453 i = ni(j)
454 ypct(j) = pctsrf(i, nsrf)
455 yts(j) = ts(i, nsrf)
456 ytslab(i) = tslab(i)
457 ysnow(j) = snow(i, nsrf)
458 yqsurf(j) = qsurf(i, nsrf)
459 yalb(j) = albe(i, nsrf)
460 yalblw(j) = alblw(i, nsrf)
461 yrain_f(j) = rain_f(i)
462 ysnow_f(j) = snow_f(i)
463 yagesno(j) = agesno(i, nsrf)
464 yfder(j) = fder(i)
465 ytaux(j) = flux_u(i, 1, nsrf)
466 ytauy(j) = flux_v(i, 1, nsrf)
467 ysolsw(j) = solsw(i, nsrf)
468 ysollw(j) = sollw(i, nsrf)
469 ysollwdown(j) = sollwdown(i)
470 yrugos(j) = rugos(i, nsrf)
471 yrugoro(j) = rugoro(i)
472 yu1(j) = u1lay(i)
473 yv1(j) = v1lay(i)
474 yrads(j) = ysolsw(j) + ysollw(j)
475 ypaprs(j, klev+1) = paprs(i, klev+1)
476 y_run_off_lic_0(j) = run_off_lic_0(i)
477 yu10mx(j) = u10m(i, nsrf)
478 yu10my(j) = v10m(i, nsrf)
479 ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
480 END DO
481
482 ! IF bucket model for continent, copy soil water content
483 IF (nsrf==is_ter .AND. .NOT. ok_veget) THEN
484 DO j = 1, knon
485 i = ni(j)
486 yqsol(j) = qsol(i)
487 END DO
488 ELSE
489 yqsol = 0.
490 END IF
491 !$$$ PB ajour pour soil
492 DO k = 1, nsoilmx
493 DO j = 1, knon
494 i = ni(j)
495 ytsoil(j, k) = ftsoil(i, k, nsrf)
496 END DO
497 END DO
498 DO k = 1, klev
499 DO j = 1, knon
500 i = ni(j)
501 ypaprs(j, k) = paprs(i, k)
502 ypplay(j, k) = pplay(i, k)
503 ydelp(j, k) = delp(i, k)
504 yu(j, k) = u(i, k)
505 yv(j, k) = v(i, k)
506 yt(j, k) = t(i, k)
507 yq(j, k) = q(i, k)
508 END DO
509 END DO
510
511 ! calculer Cdrag et les coefficients d'echange
512 CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&
513 yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)
514 !IM 081204 BEG
515 !CR test
516 IF (iflag_pbl==1) THEN
517 !IM 081204 END
518 CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
519 DO k = 1, klev
520 DO i = 1, knon
521 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
522 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
523 END DO
524 END DO
525 END IF
526
527 !IM cf JLD : on seuille ycoefm et ycoefh
528 IF (nsrf==is_oce) THEN
529 DO j = 1, knon
530 ! ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)
531 ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)
532 ! ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)
533 ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)
534 END DO
535 END IF
536
537 !IM: 261103
538 IF (ok_kzmin) THEN
539 !IM cf FH: 201103 BEG
540 ! Calcul d'une diffusion minimale pour les conditions tres stables.
541 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, &
542 ycoefm0, ycoefh0)
543
544 IF (1==1) THEN
545 DO k = 1, klev
546 DO i = 1, knon
547 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
548 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
549 END DO
550 END DO
551 END IF
552 !IM cf FH: 201103 END
553 !IM: 261103
554 END IF !ok_kzmin
555
556 IF (iflag_pbl>=3) THEN
557 ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et Frédéric Hourdin
558 yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &
559 1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg
560 DO k = 2, klev
561 yzlay(1:knon, k) = yzlay(1:knon, k-1) &
562 + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
563 / ypaprs(1:knon, k) &
564 * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
565 END DO
566 DO k = 1, klev
567 yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
568 / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
569 END DO
570 yzlev(1:knon, 1) = 0.
571 yzlev(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1:knon, klev-1)
572 DO k = 2, klev
573 yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
574 END DO
575 DO k = 1, klev + 1
576 DO j = 1, knon
577 i = ni(j)
578 yq2(j, k) = q2(i, k, nsrf)
579 END DO
580 END DO
581
582 ! Bug introduit volontairement pour converger avec les resultats
583 ! du papier sur les thermiques.
584 IF (1==1) THEN
585 y_cd_m(1:knon) = ycoefm(1:knon, 1)
586 y_cd_h(1:knon) = ycoefh(1:knon, 1)
587 ELSE
588 y_cd_h(1:knon) = ycoefm(1:knon, 1)
589 y_cd_m(1:knon) = ycoefh(1:knon, 1)
590 END IF
591 CALL ustarhb(knon, yu, yv, y_cd_m, yustar)
592
593 IF (prt_level>9) THEN
594 PRINT *, 'USTAR = ', yustar
595 END IF
596
597 ! iflag_pbl peut etre utilise comme longuer de melange
598
599 IF (iflag_pbl>=11) THEN
600 CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
601 yu, yv, yteta, y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, &
602 iflag_pbl)
603 ELSE
604 CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, &
605 yv, yteta, y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
606 END IF
607
608 ycoefm(1:knon, 1) = y_cd_m(1:knon)
609 ycoefh(1:knon, 1) = y_cd_h(1:knon)
610 ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)
611 ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)
612 END IF
613
614 ! calculer la diffusion des vitesses "u" et "v"
615 CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &
616 ydelp, y_d_u, y_flux_u)
617 CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &
618 ydelp, y_d_v, y_flux_v)
619
620 ! pour le couplage
621 ytaux = y_flux_u(:, 1)
622 ytauy = y_flux_v(:, 1)
623
624 ! FH modif sur le cdrag temperature
625 !$$$PB : déplace dans clcdrag
626 !$$$ do i=1, knon
627 !$$$ ycoefh(i, 1)=ycoefm(i, 1)*0.8
628 !$$$ enddo
629
630 ! calculer la diffusion de "q" et de "h"
631 CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&
632 cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&
633 yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&
634 yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&
635 ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &
636 yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&
637 yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&
638 yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&
639 y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&
640 ytslab, y_seaice)
641
642 ! calculer la longueur de rugosite sur ocean
643 yrugm = 0.
644 IF (nsrf==is_oce) THEN
645 DO j = 1, knon
646 yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
647 0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))
648 yrugm(j) = max(1.5E-05, yrugm(j))
649 END DO
650 END IF
651 DO j = 1, knon
652 y_dflux_t(j) = y_dflux_t(j)*ypct(j)
653 y_dflux_q(j) = y_dflux_q(j)*ypct(j)
654 yu1(j) = yu1(j)*ypct(j)
655 yv1(j) = yv1(j)*ypct(j)
656 END DO
657
658 DO k = 1, klev
659 DO j = 1, knon
660 i = ni(j)
661 ycoefh(j, k) = ycoefh(j, k)*ypct(j)
662 ycoefm(j, k) = ycoefm(j, k)*ypct(j)
663 y_d_t(j, k) = y_d_t(j, k)*ypct(j)
664 y_d_q(j, k) = y_d_q(j, k)*ypct(j)
665 !§§§ PB
666 flux_t(i, k, nsrf) = y_flux_t(j, k)
667 flux_q(i, k, nsrf) = y_flux_q(j, k)
668 flux_u(i, k, nsrf) = y_flux_u(j, k)
669 flux_v(i, k, nsrf) = y_flux_v(j, k)
670 !$$$ PB y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)
671 !$$$ PB y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)
672 y_d_u(j, k) = y_d_u(j, k)*ypct(j)
673 y_d_v(j, k) = y_d_v(j, k)*ypct(j)
674 !$$$ PB y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)
675 !$$$ PB y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)
676 END DO
677 END DO
678
679 evap(:, nsrf) = -flux_q(:, 1, nsrf)
680
681 albe(:, nsrf) = 0.
682 alblw(:, nsrf) = 0.
683 snow(:, nsrf) = 0.
684 qsurf(:, nsrf) = 0.
685 rugos(:, nsrf) = 0.
686 fluxlat(:, nsrf) = 0.
687 DO j = 1, knon
688 i = ni(j)
689 d_ts(i, nsrf) = y_d_ts(j)
690 albe(i, nsrf) = yalb(j)
691 alblw(i, nsrf) = yalblw(j)
692 snow(i, nsrf) = ysnow(j)
693 qsurf(i, nsrf) = yqsurf(j)
694 rugos(i, nsrf) = yz0_new(j)
695 fluxlat(i, nsrf) = yfluxlat(j)
696 !$$$ pb rugmer(i) = yrugm(j)
697 IF (nsrf==is_oce) THEN
698 rugmer(i) = yrugm(j)
699 rugos(i, nsrf) = yrugm(j)
700 END IF
701 !IM cf JLD ??
702 agesno(i, nsrf) = yagesno(j)
703 fqcalving(i, nsrf) = y_fqcalving(j)
704 ffonte(i, nsrf) = y_ffonte(j)
705 cdragh(i) = cdragh(i) + ycoefh(j, 1)
706 cdragm(i) = cdragm(i) + ycoefm(j, 1)
707 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
708 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
709 zu1(i) = zu1(i) + yu1(j)
710 zv1(i) = zv1(i) + yv1(j)
711 END DO
712 IF (nsrf==is_ter) THEN
713 DO j = 1, knon
714 i = ni(j)
715 qsol(i) = yqsol(j)
716 END DO
717 END IF
718 IF (nsrf==is_lic) THEN
719 DO j = 1, knon
720 i = ni(j)
721 run_off_lic_0(i) = y_run_off_lic_0(j)
722 END DO
723 END IF
724 !$$$ PB ajout pour soil
725 ftsoil(:, :, nsrf) = 0.
726 DO k = 1, nsoilmx
727 DO j = 1, knon
728 i = ni(j)
729 ftsoil(i, k, nsrf) = ytsoil(j, k)
730 END DO
731 END DO
732
733 DO j = 1, knon
734 i = ni(j)
735 DO k = 1, klev
736 d_t(i, k) = d_t(i, k) + y_d_t(j, k)
737 d_q(i, k) = d_q(i, k) + y_d_q(j, k)
738 !$$$ PB flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)
739 !$$$ flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k)
740 d_u(i, k) = d_u(i, k) + y_d_u(j, k)
741 d_v(i, k) = d_v(i, k) + y_d_v(j, k)
742 !$$$ PB flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)
743 !$$$ flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)
744 zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)
745 END DO
746 END DO
747
748 !cc diagnostic t, q a 2m et u, v a 10m
749
750 DO j = 1, knon
751 i = ni(j)
752 uzon(j) = yu(j, 1) + y_d_u(j, 1)
753 vmer(j) = yv(j, 1) + y_d_v(j, 1)
754 tair1(j) = yt(j, 1) + y_d_t(j, 1)
755 qair1(j) = yq(j, 1) + y_d_q(j, 1)
756 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
757 1)))*(ypaprs(j, 1)-ypplay(j, 1))
758 tairsol(j) = yts(j) + y_d_ts(j)
759 rugo1(j) = yrugos(j)
760 IF (nsrf==is_oce) THEN
761 rugo1(j) = rugos(i, nsrf)
762 END IF
763 psfce(j) = ypaprs(j, 1)
764 patm(j) = ypplay(j, 1)
765
766 qairsol(j) = yqsurf(j)
767 END DO
768
769 CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &
770 tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &
771 yu10m, yustar)
772 !IM 081204 END
773
774 DO j = 1, knon
775 i = ni(j)
776 t2m(i, nsrf) = yt2m(j)
777 q2m(i, nsrf) = yq2m(j)
778
779 ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
780 u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
781 v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
782
783 END DO
784
785 DO i = 1, knon
786 y_cd_h(i) = ycoefh(i, 1)
787 y_cd_m(i) = ycoefm(i, 1)
788 END DO
789 CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
790 y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
791 ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
792
793 DO j = 1, knon
794 i = ni(j)
795 pblh(i, nsrf) = ypblh(j)
796 plcl(i, nsrf) = ylcl(j)
797 capcl(i, nsrf) = ycapcl(j)
798 oliqcl(i, nsrf) = yoliqcl(j)
799 cteicl(i, nsrf) = ycteicl(j)
800 pblt(i, nsrf) = ypblt(j)
801 therm(i, nsrf) = ytherm(j)
802 trmb1(i, nsrf) = ytrmb1(j)
803 trmb2(i, nsrf) = ytrmb2(j)
804 trmb3(i, nsrf) = ytrmb3(j)
805 END DO
806
807 DO j = 1, knon
808 DO k = 1, klev + 1
809 i = ni(j)
810 q2(i, k, nsrf) = yq2(j, k)
811 END DO
812 END DO
813 !IM "slab" ocean
814 IF (nsrf==is_oce) THEN
815 DO j = 1, knon
816 ! on projette sur la grille globale
817 i = ni(j)
818 IF (pctsrf_new(i, is_oce)>epsfra) THEN
819 flux_o(i) = y_flux_o(j)
820 ELSE
821 flux_o(i) = 0.
822 END IF
823 END DO
824 END IF
825
826 IF (nsrf==is_sic) THEN
827 DO j = 1, knon
828 i = ni(j)
829 ! On pondère lorsque l'on fait le bilan au sol :
830 ! flux_g(i) = y_flux_g(j)*ypct(j)
831 IF (pctsrf_new(i, is_sic)>epsfra) THEN
832 flux_g(i) = y_flux_g(j)
833 ELSE
834 flux_g(i) = 0.
835 END IF
836 END DO
837
838 END IF
839 !nsrf.EQ.is_sic
840 IF (ocean=='slab ') THEN
841 IF (nsrf==is_oce) THEN
842 tslab(1:klon) = ytslab(1:klon)
843 seaice(1:klon) = y_seaice(1:klon)
844 !nsrf
845 END IF
846 !OCEAN
847 END IF
848 END DO
849
850 ! On utilise les nouvelles surfaces
851 ! A rajouter: conservation de l'albedo
852
853 rugos(:, is_oce) = rugmer
854 pctsrf = pctsrf_new
855
856 END SUBROUTINE clmain
857
858 end module clmain_m

  ViewVC Help
Powered by ViewVC 1.1.21