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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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

  ViewVC Help
Powered by ViewVC 1.1.21