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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 25 - (show annotations)
Fri Mar 5 16:43:45 2010 UTC (14 years, 1 month ago) by guez
Original Path: trunk/libf/phylmd/clmain.f90
File size: 27774 byte(s)
Simplified "etat0_lim.sh" and "gcm.sh" because the full versions
depended on personal arrangements for directories and machines.

Translated included files into modules. Encapsulated procedures into modules.

Moved variables from module "comgeom" to local variables of
"inigeom". Deleted some unused variables in "comgeom".

Moved variable "day_ini" from module "temps" to module "dynetat0_m".

Removed useless test on variable "time" and useless "close" statement
in procedure "leapfrog".

Removed useless call to "inigeom" in procedure "limit".

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

  ViewVC Help
Powered by ViewVC 1.1.21