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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 15 - (show annotations)
Fri Aug 1 15:24:12 2008 UTC (15 years, 8 months ago) by guez
File size: 27540 byte(s)
-- Minor modification of input/output:

Added variable "Sigma_O3_Royer" to "histday.nc". "ecrit_day" is not
modified in "physiq". Removed variables "pyu1", "pyv1", "ftsol1",
"ftsol2", "ftsol3", "ftsol4", "psrf1", "psrf2", "psrf3", "psrf4"
"mfu", "mfd", "en_u", "en_d", "de_d", "de_u", "coefh" from
"histrac.nc".

Variable "raz_date" of module "conf_gcm_m" has logical type instead of
integer type.

-- Should not change any result at run time:

Modified calls to "IOIPSL_Lionel" procedures because the interfaces of
these procedures have been simplified.

Changed name of variable in module "start_init_orog_m": "masque" to
"mask".

Created a module containing procedure "phyredem".

Removed arguments "punjours", "pdayref" and "ptimestep" of procedure
"iniphysiq".

Renamed procedure "gr_phy_write" to "gr_phy_write_2d". Created
procedure "gr_phy_write_3d".

Removed procedures "ini_undefstd", "moy_undefSTD", "calcul_STDlev",
"calcul_divers".

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

  ViewVC Help
Powered by ViewVC 1.1.21