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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 30 - (show annotations)
Thu Apr 1 09:07:28 2010 UTC (14 years ago) by guez
Original Path: trunk/libf/phylmd/clmain.f90
File size: 27812 byte(s)
Imported Source files of the external library "IOIPSL_Lionel" into
"libf/IOIPSL".

Split "cray.f90" into "scopy.f90" and "ssum.f90".

Rewrote "leapfrog" in order to have a clearer algorithmic structure.

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 histcom, ONLY : histbeg_totreg, histdef, histend, histsync
89 use histwrite_m, only: histwrite
90 use calendar, ONLY : ymds2ju
91 USE dimens_m, ONLY : iim, jjm
92 USE indicesol, ONLY : epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
93 USE dimphy, ONLY : klev, klon, zmasq
94 USE dimsoil, ONLY : nsoilmx
95 USE temps, ONLY : annee_ref, itau_phy
96 USE dynetat0_m, ONLY : day_ini
97 USE iniprint, ONLY : prt_level
98 USE yomcst, ONLY : rd, rg, rkappa
99 USE conf_phys_m, ONLY : iflag_pbl
100 USE gath_cpl, ONLY : gath2cpl
101
102 IMPLICIT NONE
103
104 REAL, INTENT (IN) :: dtime
105 REAL date0
106 INTEGER, INTENT (IN) :: itap
107 REAL t(klon, klev), q(klon, klev)
108 REAL u(klon, klev), v(klon, klev)
109 REAL, INTENT (IN) :: paprs(klon, klev+1)
110 REAL, INTENT (IN) :: pplay(klon, klev)
111 REAL, INTENT (IN) :: rlon(klon), rlat(klon)
112 REAL cufi(klon), cvfi(klon)
113 REAL d_t(klon, klev), d_q(klon, klev)
114 REAL d_u(klon, klev), d_v(klon, klev)
115 REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
116 REAL dflux_t(klon), dflux_q(klon)
117 !IM "slab" ocean
118 REAL flux_o(klon), flux_g(klon)
119 REAL y_flux_o(klon), y_flux_g(klon)
120 REAL tslab(klon), ytslab(klon)
121 REAL seaice(klon), y_seaice(klon)
122 REAL y_fqcalving(klon), y_ffonte(klon)
123 REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
124 REAL run_off_lic_0(klon), y_run_off_lic_0(klon)
125
126 REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
127 REAL rugmer(klon), agesno(klon, nbsrf)
128 REAL, INTENT (IN) :: rugoro(klon)
129 REAL cdragh(klon), cdragm(klon)
130 ! jour de l'annee en cours
131 INTEGER jour
132 REAL rmu0(klon) ! cosinus de l'angle solaire zenithal
133 ! taux CO2 atmosphere
134 REAL co2_ppm
135 LOGICAL, INTENT (IN) :: debut
136 LOGICAL, INTENT (IN) :: lafin
137 LOGICAL ok_veget
138 CHARACTER (len=*), INTENT (IN) :: ocean
139 INTEGER npas, nexca
140
141 REAL pctsrf(klon, nbsrf)
142 REAL ts(klon, nbsrf)
143 REAL d_ts(klon, nbsrf)
144 REAL snow(klon, nbsrf)
145 REAL qsurf(klon, nbsrf)
146 REAL evap(klon, nbsrf)
147 REAL albe(klon, nbsrf)
148 REAL alblw(klon, nbsrf)
149
150 REAL fluxlat(klon, nbsrf)
151
152 REAL rain_f(klon), snow_f(klon)
153 REAL fder(klon)
154
155 REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)
156 REAL rugos(klon, nbsrf)
157 ! la nouvelle repartition des surfaces sortie de l'interface
158 REAL pctsrf_new(klon, nbsrf)
159
160 REAL zcoefh(klon, klev)
161 REAL zu1(klon)
162 REAL zv1(klon)
163
164 !$$$ PB ajout pour soil
165 LOGICAL, INTENT (IN) :: soil_model
166 !IM ajout seuils cdrm, cdrh
167 REAL cdmmax, cdhmax
168
169 REAL ksta, ksta_ter
170 LOGICAL ok_kzmin
171
172 REAL ftsoil(klon, nsoilmx, nbsrf)
173 REAL ytsoil(klon, nsoilmx)
174 REAL qsol(klon)
175
176 EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
177
178 REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
179 REAL yalb(klon)
180 REAL yalblw(klon)
181 REAL yu1(klon), yv1(klon)
182 REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
183 REAL yrain_f(klon), ysnow_f(klon)
184 REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)
185 REAL yfder(klon), ytaux(klon), ytauy(klon)
186 REAL yrugm(klon), yrads(klon), yrugoro(klon)
187
188 REAL yfluxlat(klon)
189
190 REAL y_d_ts(klon)
191 REAL y_d_t(klon, klev), y_d_q(klon, klev)
192 REAL y_d_u(klon, klev), y_d_v(klon, klev)
193 REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
194 REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
195 REAL y_dflux_t(klon), y_dflux_q(klon)
196 REAL ycoefh(klon, klev), ycoefm(klon, klev)
197 REAL yu(klon, klev), yv(klon, klev)
198 REAL yt(klon, klev), yq(klon, klev)
199 REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
200
201 LOGICAL ok_nonloc
202 PARAMETER (ok_nonloc=.FALSE.)
203 REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
204
205 !IM 081204 hcl_Anne ? BEG
206 REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
207 REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
208 REAL ykmq(klon, klev+1)
209 REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)
210 REAL q2diag(klon, klev+1)
211 !IM 081204 hcl_Anne ? END
212
213 REAL u1lay(klon), v1lay(klon)
214 REAL delp(klon, klev)
215 INTEGER i, k, nsrf
216
217 INTEGER ni(klon), knon, j
218 ! Introduction d'une variable "pourcentage potentiel" pour tenir compte
219 ! des eventuelles apparitions et/ou disparitions de la glace de mer
220 REAL pctsrf_pot(klon, nbsrf)
221
222 REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
223
224 ! maf pour sorties IOISPL en cas de debugagage
225
226 CHARACTER (80) cldebug
227 SAVE cldebug
228 CHARACTER (8) cl_surf(nbsrf)
229 SAVE cl_surf
230 INTEGER nhoridbg, nidbg
231 SAVE nhoridbg, nidbg
232 INTEGER ndexbg(iim*(jjm+1))
233 REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian
234 REAL tabindx(klon)
235 REAL debugtab(iim, jjm+1)
236 LOGICAL first_appel
237 SAVE first_appel
238 DATA first_appel/ .TRUE./
239 LOGICAL :: debugindex = .FALSE.
240 INTEGER idayref
241 REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
242 REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
243
244 REAL yt2m(klon), yq2m(klon), yu10m(klon)
245 REAL yustar(klon)
246 ! -- LOOP
247 REAL yu10mx(klon)
248 REAL yu10my(klon)
249 REAL ywindsp(klon)
250 ! -- LOOP
251
252 REAL yt10m(klon), yq10m(klon)
253 !IM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds
254 ! physiq ce qui permet de sortir les grdeurs par sous surface)
255 REAL pblh(klon, nbsrf)
256 REAL plcl(klon, nbsrf)
257 REAL capcl(klon, nbsrf)
258 REAL oliqcl(klon, nbsrf)
259 REAL cteicl(klon, nbsrf)
260 REAL pblt(klon, nbsrf)
261 REAL therm(klon, nbsrf)
262 REAL trmb1(klon, nbsrf)
263 REAL trmb2(klon, nbsrf)
264 REAL trmb3(klon, nbsrf)
265 REAL ypblh(klon)
266 REAL ylcl(klon)
267 REAL ycapcl(klon)
268 REAL yoliqcl(klon)
269 REAL ycteicl(klon)
270 REAL ypblt(klon)
271 REAL ytherm(klon)
272 REAL ytrmb1(klon)
273 REAL ytrmb2(klon)
274 REAL ytrmb3(klon)
275 REAL y_cd_h(klon), y_cd_m(klon)
276 REAL uzon(klon), vmer(klon)
277 REAL tair1(klon), qair1(klon), tairsol(klon)
278 REAL psfce(klon), patm(klon)
279
280 REAL qairsol(klon), zgeo1(klon)
281 REAL rugo1(klon)
282
283 ! utiliser un jeu de fonctions simples
284 LOGICAL zxli
285 PARAMETER (zxli=.FALSE.)
286
287 REAL zt, zqs, zdelta, zcor
288 REAL t_coup
289 PARAMETER (t_coup=273.15)
290
291 CHARACTER (len=20) :: modname = 'clmain'
292 LOGICAL check
293 PARAMETER (check=.FALSE.)
294
295 !------------------------------------------------------------
296
297 ! initialisation Anne
298 ytherm = 0.
299
300 IF (check) THEN
301 PRINT *, modname, ' klon=', klon
302 END IF
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 !$$$ flux_t(i, k) = 0.0
410 !$$$ flux_q(i, k) = 0.0
411 d_u(i, k) = 0.0
412 d_v(i, k) = 0.0
413 !$$$ flux_u(i, k) = 0.0
414 !$$$ flux_v(i, k) = 0.0
415 zcoefh(i, k) = 0.0
416 END DO
417 END DO
418 !AA IF (itr.GE.1) THEN
419 !AA DO it = 1, itr
420 !AA DO k = 1, klev
421 !AA DO i = 1, klon
422 !AA d_tr(i, k, it) = 0.0
423 !AA ENDDO
424 !AA ENDDO
425 !AA ENDDO
426 !AA ENDIF
427
428
429 ! Boucler sur toutes les sous-fractions du sol:
430
431 ! Initialisation des "pourcentages potentiels". On considere ici qu'on
432 ! peut avoir potentiellementdela glace sur tout le domaine oceanique
433 ! (a affiner)
434
435 pctsrf_pot = pctsrf
436 pctsrf_pot(:, is_oce) = 1. - zmasq
437 pctsrf_pot(:, is_sic) = 1. - zmasq
438
439 DO nsrf = 1, nbsrf
440 ! chercher les indices:
441 ni = 0
442 knon = 0
443 DO i = 1, klon
444 ! pour determiner le domaine a traiter on utilise les surfaces
445 ! "potentielles"
446 IF (pctsrf_pot(i, nsrf) > epsfra) THEN
447 knon = knon + 1
448 ni(knon) = i
449 END IF
450 END DO
451
452 IF (check) THEN
453 PRINT *, 'CLMAIN, nsrf, knon =', nsrf, knon
454 END IF
455
456 ! variables pour avoir une sortie IOIPSL des INDEX
457 IF (debugindex) THEN
458 tabindx = 0.
459 DO i = 1, knon
460 tabindx(i) = real(i)
461 END DO
462 debugtab = 0.
463 ndexbg = 0
464 CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)
465 CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
466 END IF
467
468 IF (knon==0) CYCLE
469
470 DO j = 1, knon
471 i = ni(j)
472 ypct(j) = pctsrf(i, nsrf)
473 yts(j) = ts(i, nsrf)
474 ytslab(i) = tslab(i)
475 ysnow(j) = snow(i, nsrf)
476 yqsurf(j) = qsurf(i, nsrf)
477 yalb(j) = albe(i, nsrf)
478 yalblw(j) = alblw(i, nsrf)
479 yrain_f(j) = rain_f(i)
480 ysnow_f(j) = snow_f(i)
481 yagesno(j) = agesno(i, nsrf)
482 yfder(j) = fder(i)
483 ytaux(j) = flux_u(i, 1, nsrf)
484 ytauy(j) = flux_v(i, 1, nsrf)
485 ysolsw(j) = solsw(i, nsrf)
486 ysollw(j) = sollw(i, nsrf)
487 ysollwdown(j) = sollwdown(i)
488 yrugos(j) = rugos(i, nsrf)
489 yrugoro(j) = rugoro(i)
490 yu1(j) = u1lay(i)
491 yv1(j) = v1lay(i)
492 yrads(j) = ysolsw(j) + ysollw(j)
493 ypaprs(j, klev+1) = paprs(i, klev+1)
494 y_run_off_lic_0(j) = run_off_lic_0(i)
495 yu10mx(j) = u10m(i, nsrf)
496 yu10my(j) = v10m(i, nsrf)
497 ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
498 END DO
499
500 ! IF bucket model for continent, copy soil water content
501 IF (nsrf==is_ter .AND. .NOT. ok_veget) THEN
502 DO j = 1, knon
503 i = ni(j)
504 yqsol(j) = qsol(i)
505 END DO
506 ELSE
507 yqsol = 0.
508 END IF
509 !$$$ PB ajour pour soil
510 DO k = 1, nsoilmx
511 DO j = 1, knon
512 i = ni(j)
513 ytsoil(j, k) = ftsoil(i, k, nsrf)
514 END DO
515 END DO
516 DO k = 1, klev
517 DO j = 1, knon
518 i = ni(j)
519 ypaprs(j, k) = paprs(i, k)
520 ypplay(j, k) = pplay(i, k)
521 ydelp(j, k) = delp(i, k)
522 yu(j, k) = u(i, k)
523 yv(j, k) = v(i, k)
524 yt(j, k) = t(i, k)
525 yq(j, k) = q(i, k)
526 END DO
527 END DO
528
529 ! calculer Cdrag et les coefficients d'echange
530 CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts,&
531 yrugos, yu, yv, yt, yq, yqsurf, ycoefm, ycoefh)
532 !IM 081204 BEG
533 !CR test
534 IF (iflag_pbl==1) THEN
535 !IM 081204 END
536 CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
537 DO k = 1, klev
538 DO i = 1, knon
539 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
540 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
541 END DO
542 END DO
543 END IF
544
545 !IM cf JLD : on seuille ycoefm et ycoefh
546 IF (nsrf==is_oce) THEN
547 DO j = 1, knon
548 ! ycoefm(j, 1)=min(ycoefm(j, 1), 1.1E-3)
549 ycoefm(j, 1) = min(ycoefm(j, 1), cdmmax)
550 ! ycoefh(j, 1)=min(ycoefh(j, 1), 1.1E-3)
551 ycoefh(j, 1) = min(ycoefh(j, 1), cdhmax)
552 END DO
553 END IF
554
555
556 !IM: 261103
557 IF (ok_kzmin) THEN
558 !IM cf FH: 201103 BEG
559 ! Calcul d'une diffusion minimale pour les conditions tres stables.
560 CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, ycoefm, ycoefm0, &
561 ycoefh0)
562 ! call dump2d(iim, jjm-1, ycoefm(2:klon-1, 2), 'KZ ')
563 ! call dump2d(iim, jjm-1, ycoefm0(2:klon-1, 2), 'KZMIN ')
564
565 IF (1==1) THEN
566 DO k = 1, klev
567 DO i = 1, knon
568 ycoefm(i, k) = max(ycoefm(i, k), ycoefm0(i, k))
569 ycoefh(i, k) = max(ycoefh(i, k), ycoefh0(i, k))
570 END DO
571 END DO
572 END IF
573 !IM cf FH: 201103 END
574 !IM: 261103
575 END IF !ok_kzmin
576
577 IF (iflag_pbl>=3) THEN
578
579 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
580 ! MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
581 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
582
583 yzlay(1:knon, 1) = rd*yt(1:knon, 1)/(0.5*(ypaprs(1:knon, &
584 1)+ypplay(1:knon, 1)))*(ypaprs(1:knon, 1)-ypplay(1:knon, 1))/rg
585 DO k = 2, klev
586 yzlay(1:knon, k) = yzlay(1:knon, k-1) &
587 + rd*0.5*(yt(1:knon, k-1) +yt(1: knon, k)) &
588 / ypaprs(1:knon, k) *(ypplay(1:knon, k-1)-ypplay(1:knon, k))/ &
589 rg
590 END DO
591 DO k = 1, klev
592 yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
593 / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
594 END DO
595 yzlev(1:knon, 1) = 0.
596 yzlev(1:knon, klev+1) = 2.*yzlay(1:knon, klev) - yzlay(1:knon, klev-1)
597 DO k = 2, klev
598 yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
599 END DO
600 DO k = 1, klev + 1
601 DO j = 1, knon
602 i = ni(j)
603 yq2(j, k) = q2(i, k, nsrf)
604 END DO
605 END DO
606
607
608 ! Bug introduit volontairement pour converger avec les resultats
609 ! du papier sur les thermiques.
610 IF (1==1) THEN
611 y_cd_m(1:knon) = ycoefm(1:knon, 1)
612 y_cd_h(1:knon) = ycoefh(1:knon, 1)
613 ELSE
614 y_cd_h(1:knon) = ycoefm(1:knon, 1)
615 y_cd_m(1:knon) = ycoefh(1:knon, 1)
616 END IF
617 CALL ustarhb(knon, yu, yv, y_cd_m, yustar)
618
619 IF (prt_level>9) THEN
620 PRINT *, 'USTAR = ', yustar
621 END IF
622
623 ! iflag_pbl peut etre utilise comme longuer de melange
624
625 IF (iflag_pbl>=11) THEN
626 CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, yv, yteta, &
627 y_cd_m, yq2, q2diag, ykmm, ykmn, yustar, iflag_pbl)
628 ELSE
629 CALL yamada4(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, yu, yv, yteta, &
630 y_cd_m, yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
631 END IF
632
633 ycoefm(1:knon, 1) = y_cd_m(1:knon)
634 ycoefh(1:knon, 1) = y_cd_h(1:knon)
635 ycoefm(1:knon, 2:klev) = ykmm(1:knon, 2:klev)
636 ycoefh(1:knon, 2:klev) = ykmn(1:knon, 2:klev)
637
638
639 END IF
640
641 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
642 ! calculer la diffusion des vitesses "u" et "v"
643 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
644
645 CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yu, ypaprs, ypplay, &
646 ydelp, y_d_u, y_flux_u)
647 CALL clvent(knon, dtime, yu1, yv1, ycoefm, yt, yv, ypaprs, ypplay, &
648 ydelp, y_d_v, y_flux_v)
649
650 ! pour le couplage
651 ytaux = y_flux_u(:, 1)
652 ytauy = y_flux_v(:, 1)
653
654 ! FH modif sur le cdrag temperature
655 !$$$PB : déplace dans clcdrag
656 !$$$ do i=1, knon
657 !$$$ ycoefh(i, 1)=ycoefm(i, 1)*0.8
658 !$$$ enddo
659
660 ! calculer la diffusion de "q" et de "h"
661 CALL clqh(dtime, itap, date0, jour, debut, lafin, rlon, rlat,&
662 cufi, cvfi, knon, nsrf, ni, pctsrf, soil_model, ytsoil,&
663 yqsol, ok_veget, ocean, npas, nexca, rmu0, co2_ppm, yrugos,&
664 yrugoro, yu1, yv1, ycoefh, yt, yq, yts, ypaprs, ypplay,&
665 ydelp, yrads, yalb, yalblw, ysnow, yqsurf, yrain_f, ysnow_f, &
666 yfder, ytaux, ytauy, ywindsp, ysollw, ysollwdown, ysolsw,&
667 yfluxlat, pctsrf_new, yagesno, y_d_t, y_d_q, y_d_ts,&
668 yz0_new, y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,&
669 y_fqcalving, y_ffonte, y_run_off_lic_0, y_flux_o, y_flux_g,&
670 ytslab, y_seaice)
671
672 ! calculer la longueur de rugosite sur ocean
673 yrugm = 0.
674 IF (nsrf==is_oce) THEN
675 DO j = 1, knon
676 yrugm(j) = 0.018*ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
677 0.11*14E-6/sqrt(ycoefm(j, 1)*(yu1(j)**2+yv1(j)**2))
678 yrugm(j) = max(1.5E-05, yrugm(j))
679 END DO
680 END IF
681 DO j = 1, knon
682 y_dflux_t(j) = y_dflux_t(j)*ypct(j)
683 y_dflux_q(j) = y_dflux_q(j)*ypct(j)
684 yu1(j) = yu1(j)*ypct(j)
685 yv1(j) = yv1(j)*ypct(j)
686 END DO
687
688 DO k = 1, klev
689 DO j = 1, knon
690 i = ni(j)
691 ycoefh(j, k) = ycoefh(j, k)*ypct(j)
692 ycoefm(j, k) = ycoefm(j, k)*ypct(j)
693 y_d_t(j, k) = y_d_t(j, k)*ypct(j)
694 y_d_q(j, k) = y_d_q(j, k)*ypct(j)
695 !§§§ PB
696 flux_t(i, k, nsrf) = y_flux_t(j, k)
697 flux_q(i, k, nsrf) = y_flux_q(j, k)
698 flux_u(i, k, nsrf) = y_flux_u(j, k)
699 flux_v(i, k, nsrf) = y_flux_v(j, k)
700 !$$$ PB y_flux_t(j, k) = y_flux_t(j, k) * ypct(j)
701 !$$$ PB y_flux_q(j, k) = y_flux_q(j, k) * ypct(j)
702 y_d_u(j, k) = y_d_u(j, k)*ypct(j)
703 y_d_v(j, k) = y_d_v(j, k)*ypct(j)
704 !$$$ PB y_flux_u(j, k) = y_flux_u(j, k) * ypct(j)
705 !$$$ PB y_flux_v(j, k) = y_flux_v(j, k) * ypct(j)
706 END DO
707 END DO
708
709
710 evap(:, nsrf) = -flux_q(:, 1, nsrf)
711
712 albe(:, nsrf) = 0.
713 alblw(:, nsrf) = 0.
714 snow(:, nsrf) = 0.
715 qsurf(:, nsrf) = 0.
716 rugos(:, nsrf) = 0.
717 fluxlat(:, nsrf) = 0.
718 DO j = 1, knon
719 i = ni(j)
720 d_ts(i, nsrf) = y_d_ts(j)
721 albe(i, nsrf) = yalb(j)
722 alblw(i, nsrf) = yalblw(j)
723 snow(i, nsrf) = ysnow(j)
724 qsurf(i, nsrf) = yqsurf(j)
725 rugos(i, nsrf) = yz0_new(j)
726 fluxlat(i, nsrf) = yfluxlat(j)
727 !$$$ pb rugmer(i) = yrugm(j)
728 IF (nsrf==is_oce) THEN
729 rugmer(i) = yrugm(j)
730 rugos(i, nsrf) = yrugm(j)
731 END IF
732 !IM cf JLD ??
733 agesno(i, nsrf) = yagesno(j)
734 fqcalving(i, nsrf) = y_fqcalving(j)
735 ffonte(i, nsrf) = y_ffonte(j)
736 cdragh(i) = cdragh(i) + ycoefh(j, 1)
737 cdragm(i) = cdragm(i) + ycoefm(j, 1)
738 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
739 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
740 zu1(i) = zu1(i) + yu1(j)
741 zv1(i) = zv1(i) + yv1(j)
742 END DO
743 IF (nsrf==is_ter) THEN
744 DO j = 1, knon
745 i = ni(j)
746 qsol(i) = yqsol(j)
747 END DO
748 END IF
749 IF (nsrf==is_lic) THEN
750 DO j = 1, knon
751 i = ni(j)
752 run_off_lic_0(i) = y_run_off_lic_0(j)
753 END DO
754 END IF
755 !$$$ PB ajout pour soil
756 ftsoil(:, :, nsrf) = 0.
757 DO k = 1, nsoilmx
758 DO j = 1, knon
759 i = ni(j)
760 ftsoil(i, k, nsrf) = ytsoil(j, k)
761 END DO
762 END DO
763
764 DO j = 1, knon
765 i = ni(j)
766 DO k = 1, klev
767 d_t(i, k) = d_t(i, k) + y_d_t(j, k)
768 d_q(i, k) = d_q(i, k) + y_d_q(j, k)
769 !$$$ PB flux_t(i, k) = flux_t(i, k) + y_flux_t(j, k)
770 !$$$ flux_q(i, k) = flux_q(i, k) + y_flux_q(j, k)
771 d_u(i, k) = d_u(i, k) + y_d_u(j, k)
772 d_v(i, k) = d_v(i, k) + y_d_v(j, k)
773 !$$$ PB flux_u(i, k) = flux_u(i, k) + y_flux_u(j, k)
774 !$$$ flux_v(i, k) = flux_v(i, k) + y_flux_v(j, k)
775 zcoefh(i, k) = zcoefh(i, k) + ycoefh(j, k)
776 END DO
777 END DO
778
779
780 !cc diagnostic t, q a 2m et u, v a 10m
781
782 DO j = 1, knon
783 i = ni(j)
784 uzon(j) = yu(j, 1) + y_d_u(j, 1)
785 vmer(j) = yv(j, 1) + y_d_v(j, 1)
786 tair1(j) = yt(j, 1) + y_d_t(j, 1)
787 qair1(j) = yq(j, 1) + y_d_q(j, 1)
788 zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
789 1)))*(ypaprs(j, 1)-ypplay(j, 1))
790 tairsol(j) = yts(j) + y_d_ts(j)
791 rugo1(j) = yrugos(j)
792 IF (nsrf==is_oce) THEN
793 rugo1(j) = rugos(i, nsrf)
794 END IF
795 psfce(j) = ypaprs(j, 1)
796 patm(j) = ypplay(j, 1)
797
798 qairsol(j) = yqsurf(j)
799 END DO
800
801 CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, zgeo1, &
802 tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, yt10m, yq10m, &
803 yu10m, yustar)
804 !IM 081204 END
805
806 DO j = 1, knon
807 i = ni(j)
808 t2m(i, nsrf) = yt2m(j)
809 q2m(i, nsrf) = yq2m(j)
810
811 ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
812 u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
813 v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
814
815 END DO
816
817 !IM cf AM : pbl, HBTM
818 DO i = 1, knon
819 y_cd_h(i) = ycoefh(i, 1)
820 y_cd_m(i) = ycoefm(i, 1)
821 END DO
822 ! print*, 'appel hbtm2'
823 CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, y_flux_t, &
824 y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, ycteicl, ypblt, ytherm, &
825 ytrmb1, ytrmb2, ytrmb3, ylcl)
826 ! print*, 'fin hbtm2'
827
828 DO j = 1, knon
829 i = ni(j)
830 pblh(i, nsrf) = ypblh(j)
831 plcl(i, nsrf) = ylcl(j)
832 capcl(i, nsrf) = ycapcl(j)
833 oliqcl(i, nsrf) = yoliqcl(j)
834 cteicl(i, nsrf) = ycteicl(j)
835 pblt(i, nsrf) = ypblt(j)
836 therm(i, nsrf) = ytherm(j)
837 trmb1(i, nsrf) = ytrmb1(j)
838 trmb2(i, nsrf) = ytrmb2(j)
839 trmb3(i, nsrf) = ytrmb3(j)
840 END DO
841
842
843 DO j = 1, knon
844 DO k = 1, klev + 1
845 i = ni(j)
846 q2(i, k, nsrf) = yq2(j, k)
847 END DO
848 END DO
849 !IM "slab" ocean
850 IF (nsrf==is_oce) THEN
851 DO j = 1, knon
852 ! on projette sur la grille globale
853 i = ni(j)
854 IF (pctsrf_new(i, is_oce)>epsfra) THEN
855 flux_o(i) = y_flux_o(j)
856 ELSE
857 flux_o(i) = 0.
858 END IF
859 END DO
860 END IF
861
862 IF (nsrf==is_sic) THEN
863 DO j = 1, knon
864 i = ni(j)
865 !IM 230604 on pondere lorsque l'on fait le bilan au sol : flux_g(i) = y_flux_g(j)*ypct(j)
866 IF (pctsrf_new(i, is_sic)>epsfra) THEN
867 flux_g(i) = y_flux_g(j)
868 ELSE
869 flux_g(i) = 0.
870 END IF
871 END DO
872
873 END IF
874 !nsrf.EQ.is_sic
875 IF (ocean=='slab ') THEN
876 IF (nsrf==is_oce) THEN
877 tslab(1:klon) = ytslab(1:klon)
878 seaice(1:klon) = y_seaice(1:klon)
879 !nsrf
880 END IF
881 !OCEAN
882 END IF
883 END DO
884
885 ! On utilise les nouvelles surfaces
886 ! A rajouter: conservation de l'albedo
887
888 rugos(:, is_oce) = rugmer
889 pctsrf = pctsrf_new
890
891 END SUBROUTINE clmain

  ViewVC Help
Powered by ViewVC 1.1.21