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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 98 - (hide annotations)
Tue May 13 17:23:16 2014 UTC (9 years, 10 months ago) by guez
Original Path: trunk/phylmd/clmain.f
File size: 25251 byte(s)
Split inter_barxy.f : one procedure per module, one module per
file. Grouped the files into a directory.

Split orbite.f.

Value of raz_date read from the namelist is taken into account
(resetting the step counter) even if annee_ref == anneeref and day_ref
== dayref. raz_date is no longer modified by gcm main unit. (Following
LMDZ.)

Removed argument klon of interfsur_lim. Renamed arguments lmt_alb,
lmt_rug to alb_new, z0_new (same name as corresponding actual
arguments in interfsurf_hq).

Removed argument klon of interfsurf_hq.

Removed arguments qs and d_qs of diagetpq. Were always
zero. Downgraded arguments d_qw, d_ql of diagetpq to local variables,
they were not used in physiq. Removed all computations for solid water
in diagetpq, was just zero.


Downgraded arguments fs_bound, fq_bound of diagphy to local variables,
they were not used in physiq. Encapsulated in a test on iprt all
computations in diagphy.

Removed parameter nbtr of module dimphy. Replaced it everywhere in the
program by nqmx - 2.

Removed parameter rnpb of procedure physiq. Kept the true case in
physiq and phytrac. Could not work with false case anyway.

Removed arguments klon, llm, airephy of qcheck. Removed argument ftsol
of initrrnpb, was not used.

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 72 SUBROUTINE clmain(dtime, itap, pctsrf, pctsrf_new, t, q, u, v, &
8     jour, rmu0, co2_ppm, ok_veget, ocean, ts, &
9 guez 62 soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &
10     qsol, paprs, pplay, snow, qsurf, evap, albe, alblw, fluxlat, &
11 guez 72 rain_fall, snow_f, solsw, sollw, fder, rlon, rlat, &
12     rugos, debut, agesno, rugoro, d_t, d_q, d_u, d_v, &
13 guez 62 d_ts, flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, &
14 guez 70 dflux_t, dflux_q, ycoefh, zu1, zv1, t2m, q2m, u10m, v10m, pblh, &
15 guez 62 capcl, oliqcl, cteicl, pblt, therm, trmb1, trmb2, trmb3, plcl, &
16 guez 38 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 62 ! Author: Z. X. Li (LMD/CNRS), date: 1993/08/18
20     ! Objet : interface de couche limite (diffusion verticale)
21 guez 3
22 guez 62 ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
23     ! de la couche limite pour les traceurs se fait avec "cltrac" et
24     ! ne tient pas compte de la différentiation des sous-fractions de
25     ! sol.
26 guez 3
27 guez 38 ! Pour pouvoir extraire les coefficients d'échanges et le vent
28 guez 70 ! dans la première couche, trois champs ont été créés : "ycoefh",
29 guez 62 ! "zu1" et "zv1". Nous avons moyenné les valeurs de ces trois
30     ! champs sur les quatre sous-surfaces du modèle.
31 guez 3
32 guez 49 use clqh_m, only: clqh
33 guez 62 use clvent_m, only: clvent
34 guez 47 use coefkz_m, only: coefkz
35     use coefkzmin_m, only: coefkzmin
36 guez 62 USE conf_gcm_m, ONLY: prt_level
37     USE conf_phys_m, ONLY: iflag_pbl
38     USE dimens_m, ONLY: iim, jjm
39     USE dimphy, ONLY: klev, klon, zmasq
40     USE dimsoil, ONLY: nsoilmx
41     USE dynetat0_m, ONLY: day_ini
42     USE gath_cpl, ONLY: gath2cpl
43 guez 47 use hbtm_m, only: hbtm
44 guez 62 USE histbeg_totreg_m, ONLY: histbeg_totreg
45     USE histdef_m, ONLY: histdef
46     USE histend_m, ONLY: histend
47     USE histsync_m, ONLY: histsync
48 guez 47 use histwrite_m, only: histwrite
49 guez 62 USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
50     USE suphec_m, ONLY: rd, rg, rkappa
51     USE temps, ONLY: annee_ref, itau_phy
52     use ustarhb_m, only: ustarhb
53     use vdif_kcay_m, only: vdif_kcay
54 guez 47 use yamada4_m, only: yamada4
55 guez 92 use ymds2ju_m, ONLY: ymds2ju
56 guez 15
57 guez 49 ! Arguments:
58    
59 guez 62 REAL, INTENT(IN):: dtime ! interval du temps (secondes)
60     INTEGER, INTENT(IN):: itap ! numero du pas de temps
61     REAL, INTENT(inout):: pctsrf(klon, nbsrf)
62    
63     ! la nouvelle repartition des surfaces sortie de l'interface
64     REAL, INTENT(out):: pctsrf_new(klon, nbsrf)
65    
66     REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
67     REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg/kg)
68     REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
69     INTEGER, INTENT(IN):: jour ! jour de l'annee en cours
70     REAL, intent(in):: rmu0(klon) ! cosinus de l'angle solaire zenithal
71 guez 70 REAL co2_ppm ! taux CO2 atmosphere
72     LOGICAL ok_veget
73     CHARACTER(len=*), INTENT(IN):: ocean
74     REAL ts(klon, nbsrf) ! input-R- temperature du sol (en Kelvin)
75     LOGICAL, INTENT(IN):: soil_model
76 guez 71 REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
77 guez 70 REAL ksta, ksta_ter
78     LOGICAL ok_kzmin
79     REAL ftsoil(klon, nsoilmx, nbsrf)
80     REAL qsol(klon)
81 guez 62 REAL, INTENT(IN):: paprs(klon, klev+1) ! pression a intercouche (Pa)
82     REAL, INTENT(IN):: pplay(klon, klev) ! pression au milieu de couche (Pa)
83 guez 70 REAL snow(klon, nbsrf)
84     REAL qsurf(klon, nbsrf)
85     REAL evap(klon, nbsrf)
86     REAL albe(klon, nbsrf)
87     REAL alblw(klon, nbsrf)
88    
89     REAL fluxlat(klon, nbsrf)
90    
91     REAL, intent(in):: rain_fall(klon), snow_f(klon)
92 guez 72 REAL, INTENT(IN):: solsw(klon, nbsrf), sollw(klon, nbsrf)
93 guez 70 REAL fder(klon)
94 guez 62 REAL, INTENT(IN):: rlon(klon)
95     REAL, INTENT(IN):: rlat(klon) ! latitude en degrés
96 guez 70
97     REAL rugos(klon, nbsrf)
98     ! rugos----input-R- longeur de rugosite (en m)
99    
100     LOGICAL, INTENT(IN):: debut
101     real agesno(klon, nbsrf)
102     REAL, INTENT(IN):: rugoro(klon)
103    
104 guez 38 REAL d_t(klon, klev), d_q(klon, klev)
105 guez 49 ! d_t------output-R- le changement pour "t"
106     ! d_q------output-R- le changement pour "q"
107 guez 62
108     REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
109     ! changement pour "u" et "v"
110    
111 guez 70 REAL d_ts(klon, nbsrf)
112     ! d_ts-----output-R- le changement pour "ts"
113    
114 guez 38 REAL flux_t(klon, klev, nbsrf), flux_q(klon, klev, nbsrf)
115 guez 49 ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
116     ! (orientation positive vers le bas)
117     ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
118 guez 70
119     REAL flux_u(klon, klev, nbsrf), flux_v(klon, klev, nbsrf)
120     ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
121     ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
122    
123     REAL, INTENT(out):: cdragh(klon), cdragm(klon)
124     real q2(klon, klev+1, nbsrf)
125    
126 guez 38 REAL dflux_t(klon), dflux_q(klon)
127 guez 49 ! dflux_t derive du flux sensible
128     ! dflux_q derive du flux latent
129 guez 38 !IM "slab" ocean
130 guez 70
131     REAL, intent(out):: ycoefh(klon, klev)
132     REAL, intent(out):: zu1(klon)
133     REAL zv1(klon)
134     REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
135     REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
136    
137     !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
138     ! physiq ce qui permet de sortir les grdeurs par sous surface)
139     REAL pblh(klon, nbsrf)
140     ! pblh------- HCL
141     REAL capcl(klon, nbsrf)
142     REAL oliqcl(klon, nbsrf)
143     REAL cteicl(klon, nbsrf)
144     REAL pblt(klon, nbsrf)
145     ! pblT------- T au nveau HCL
146     REAL therm(klon, nbsrf)
147     REAL trmb1(klon, nbsrf)
148     ! trmb1-------deep_cape
149     REAL trmb2(klon, nbsrf)
150     ! trmb2--------inhibition
151     REAL trmb3(klon, nbsrf)
152     ! trmb3-------Point Omega
153     REAL plcl(klon, nbsrf)
154     REAL fqcalving(klon, nbsrf), ffonte(klon, nbsrf)
155     ! ffonte----Flux thermique utilise pour fondre la neige
156     ! fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la
157     ! hauteur de neige, en kg/m2/s
158     REAL run_off_lic_0(klon)
159    
160 guez 38 REAL flux_o(klon), flux_g(klon)
161 guez 49 !IM "slab" ocean
162     ! flux_g---output-R- flux glace (pour OCEAN='slab ')
163     ! flux_o---output-R- flux ocean (pour OCEAN='slab ')
164 guez 70
165     REAL tslab(klon)
166 guez 49 ! tslab-in/output-R temperature du slab ocean (en Kelvin)
167     ! uniqmnt pour slab
168 guez 70
169     REAL seaice(klon)
170 guez 49 ! seaice---output-R- glace de mer (kg/m2) (pour OCEAN='slab ')
171 guez 15
172 guez 70 ! Local:
173 guez 15
174 guez 70 REAL y_flux_o(klon), y_flux_g(klon)
175     real ytslab(klon)
176     real y_seaice(klon)
177     REAL y_fqcalving(klon), y_ffonte(klon)
178     real y_run_off_lic_0(klon)
179 guez 15
180 guez 70 REAL rugmer(klon)
181 guez 15
182 guez 38 REAL ytsoil(klon, nsoilmx)
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 guez 49 ! on rajoute en output yu1 et yv1 qui sont les vents dans
189     ! la premiere couche
190 guez 38 REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
191     REAL yrain_f(klon), ysnow_f(klon)
192 guez 72 REAL ysollw(klon), ysolsw(klon)
193 guez 38 REAL yfder(klon), ytaux(klon), ytauy(klon)
194     REAL yrugm(klon), yrads(klon), yrugoro(klon)
195 guez 15
196 guez 38 REAL yfluxlat(klon)
197 guez 15
198 guez 38 REAL y_d_ts(klon)
199     REAL y_d_t(klon, klev), y_d_q(klon, klev)
200     REAL y_d_u(klon, klev), y_d_v(klon, klev)
201     REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
202     REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
203     REAL y_dflux_t(klon), y_dflux_q(klon)
204 guez 62 REAL coefh(klon, klev), coefm(klon, klev)
205 guez 38 REAL yu(klon, klev), yv(klon, klev)
206     REAL yt(klon, klev), yq(klon, klev)
207     REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
208 guez 15
209 guez 38 REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
210 guez 15
211 guez 38 REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
212     REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
213     REAL ykmq(klon, klev+1)
214 guez 70 REAL yq2(klon, klev+1)
215 guez 38 REAL q2diag(klon, klev+1)
216 guez 15
217 guez 38 REAL u1lay(klon), v1lay(klon)
218     REAL delp(klon, klev)
219     INTEGER i, k, nsrf
220 guez 15
221 guez 38 INTEGER ni(klon), knon, j
222 guez 40
223 guez 38 REAL pctsrf_pot(klon, nbsrf)
224 guez 40 ! "pourcentage potentiel" pour tenir compte des éventuelles
225     ! apparitions ou disparitions de la glace de mer
226 guez 15
227 guez 38 REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
228 guez 15
229 guez 38 ! maf pour sorties IOISPL en cas de debugagage
230 guez 15
231 guez 62 CHARACTER(80) cldebug
232 guez 38 SAVE cldebug
233 guez 62 CHARACTER(8) cl_surf(nbsrf)
234 guez 38 SAVE cl_surf
235     INTEGER nhoridbg, nidbg
236     SAVE nhoridbg, nidbg
237     INTEGER ndexbg(iim*(jjm+1))
238     REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian
239     REAL tabindx(klon)
240     REAL debugtab(iim, jjm+1)
241     LOGICAL first_appel
242     SAVE first_appel
243     DATA first_appel/ .TRUE./
244 guez 62 LOGICAL:: debugindex = .FALSE.
245 guez 38 INTEGER idayref
246 guez 15
247 guez 38 REAL yt2m(klon), yq2m(klon), yu10m(klon)
248     REAL yustar(klon)
249     ! -- LOOP
250     REAL yu10mx(klon)
251     REAL yu10my(klon)
252     REAL ywindsp(klon)
253     ! -- LOOP
254 guez 15
255 guez 38 REAL yt10m(klon), yq10m(klon)
256     REAL ypblh(klon)
257     REAL ylcl(klon)
258     REAL ycapcl(klon)
259     REAL yoliqcl(klon)
260     REAL ycteicl(klon)
261     REAL ypblt(klon)
262     REAL ytherm(klon)
263     REAL ytrmb1(klon)
264     REAL ytrmb2(klon)
265     REAL ytrmb3(klon)
266     REAL uzon(klon), vmer(klon)
267     REAL tair1(klon), qair1(klon), tairsol(klon)
268     REAL psfce(klon), patm(klon)
269 guez 15
270 guez 38 REAL qairsol(klon), zgeo1(klon)
271     REAL rugo1(klon)
272 guez 15
273 guez 38 ! utiliser un jeu de fonctions simples
274     LOGICAL zxli
275     PARAMETER (zxli=.FALSE.)
276 guez 15
277 guez 38 !------------------------------------------------------------
278 guez 15
279 guez 38 ytherm = 0.
280 guez 15
281 guez 38 IF (debugindex .AND. first_appel) THEN
282     first_appel = .FALSE.
283 guez 15
284 guez 38 ! initialisation sorties netcdf
285 guez 15
286 guez 38 idayref = day_ini
287 guez 40 CALL ymds2ju(annee_ref, 1, idayref, 0., zjulian)
288 guez 38 CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)
289     DO i = 1, iim
290     zx_lon(i, 1) = rlon(i+1)
291     zx_lon(i, jjm+1) = rlon(i+1)
292     END DO
293     CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlat, zx_lat)
294     cldebug = 'sous_index'
295     CALL histbeg_totreg(cldebug, zx_lon(:, 1), zx_lat(1, :), 1, &
296     iim, 1, jjm+1, itau_phy, zjulian, dtime, nhoridbg, nidbg)
297     ! no vertical axis
298     cl_surf(1) = 'ter'
299     cl_surf(2) = 'lic'
300     cl_surf(3) = 'oce'
301     cl_surf(4) = 'sic'
302     DO nsrf = 1, nbsrf
303     CALL histdef(nidbg, cl_surf(nsrf), cl_surf(nsrf), '-', iim, jjm+1, &
304     nhoridbg, 1, 1, 1, -99, 'inst', dtime, dtime)
305     END DO
306     CALL histend(nidbg)
307     CALL histsync(nidbg)
308     END IF
309 guez 15
310 guez 38 DO k = 1, klev ! epaisseur de couche
311     DO i = 1, klon
312     delp(i, k) = paprs(i, k) - paprs(i, k+1)
313     END DO
314     END DO
315     DO i = 1, klon ! vent de la premiere couche
316     zx_alf1 = 1.0
317     zx_alf2 = 1.0 - zx_alf1
318     u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2
319     v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
320     END DO
321 guez 15
322 guez 40 ! Initialization:
323     rugmer = 0.
324     cdragh = 0.
325     cdragm = 0.
326     dflux_t = 0.
327     dflux_q = 0.
328     zu1 = 0.
329     zv1 = 0.
330     ypct = 0.
331     yts = 0.
332     ysnow = 0.
333     yqsurf = 0.
334     yalb = 0.
335     yalblw = 0.
336     yrain_f = 0.
337     ysnow_f = 0.
338     yfder = 0.
339     ytaux = 0.
340     ytauy = 0.
341     ysolsw = 0.
342     ysollw = 0.
343     yrugos = 0.
344     yu1 = 0.
345     yv1 = 0.
346     yrads = 0.
347     ypaprs = 0.
348     ypplay = 0.
349     ydelp = 0.
350     yu = 0.
351     yv = 0.
352     yt = 0.
353     yq = 0.
354     pctsrf_new = 0.
355     y_flux_u = 0.
356     y_flux_v = 0.
357 guez 38 !$$ PB
358 guez 40 y_dflux_t = 0.
359     y_dflux_q = 0.
360 guez 38 ytsoil = 999999.
361     yrugoro = 0.
362     ! -- LOOP
363 guez 40 yu10mx = 0.
364     yu10my = 0.
365     ywindsp = 0.
366 guez 38 ! -- LOOP
367 guez 40 d_ts = 0.
368 guez 38 !§§§ PB
369     yfluxlat = 0.
370     flux_t = 0.
371     flux_q = 0.
372     flux_u = 0.
373     flux_v = 0.
374 guez 40 d_t = 0.
375     d_q = 0.
376     d_u = 0.
377     d_v = 0.
378 guez 70 ycoefh = 0.
379 guez 15
380 guez 38 ! Boucler sur toutes les sous-fractions du sol:
381 guez 15
382 guez 38 ! Initialisation des "pourcentages potentiels". On considère ici qu'on
383     ! peut avoir potentiellement de la glace sur tout le domaine océanique
384     ! (à affiner)
385 guez 15
386 guez 38 pctsrf_pot = pctsrf
387     pctsrf_pot(:, is_oce) = 1. - zmasq
388     pctsrf_pot(:, is_sic) = 1. - zmasq
389 guez 15
390 guez 49 loop_surface: DO nsrf = 1, nbsrf
391     ! Chercher les indices :
392 guez 38 ni = 0
393     knon = 0
394     DO i = 1, klon
395 guez 40 ! Pour déterminer le domaine à traiter, on utilise les surfaces
396 guez 38 ! "potentielles"
397     IF (pctsrf_pot(i, nsrf) > epsfra) THEN
398     knon = knon + 1
399     ni(knon) = i
400     END IF
401     END DO
402 guez 15
403 guez 38 ! variables pour avoir une sortie IOIPSL des INDEX
404     IF (debugindex) THEN
405     tabindx = 0.
406     DO i = 1, knon
407     tabindx(i) = real(i)
408     END DO
409     debugtab = 0.
410     ndexbg = 0
411     CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)
412     CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
413     END IF
414 guez 15
415 guez 62 if_knon: IF (knon /= 0) then
416 guez 38 DO j = 1, knon
417     i = ni(j)
418 guez 62 ypct(j) = pctsrf(i, nsrf)
419     yts(j) = ts(i, nsrf)
420     ytslab(i) = tslab(i)
421     ysnow(j) = snow(i, nsrf)
422     yqsurf(j) = qsurf(i, nsrf)
423     yalb(j) = albe(i, nsrf)
424     yalblw(j) = alblw(i, nsrf)
425     yrain_f(j) = rain_fall(i)
426     ysnow_f(j) = snow_f(i)
427     yagesno(j) = agesno(i, nsrf)
428     yfder(j) = fder(i)
429     ytaux(j) = flux_u(i, 1, nsrf)
430     ytauy(j) = flux_v(i, 1, nsrf)
431     ysolsw(j) = solsw(i, nsrf)
432     ysollw(j) = sollw(i, nsrf)
433     yrugos(j) = rugos(i, nsrf)
434     yrugoro(j) = rugoro(i)
435     yu1(j) = u1lay(i)
436     yv1(j) = v1lay(i)
437     yrads(j) = ysolsw(j) + ysollw(j)
438     ypaprs(j, klev+1) = paprs(i, klev+1)
439     y_run_off_lic_0(j) = run_off_lic_0(i)
440     yu10mx(j) = u10m(i, nsrf)
441     yu10my(j) = v10m(i, nsrf)
442     ywindsp(j) = sqrt(yu10mx(j)*yu10mx(j)+yu10my(j)*yu10my(j))
443 guez 38 END DO
444 guez 3
445 guez 62 ! IF bucket model for continent, copy soil water content
446     IF (nsrf == is_ter .AND. .NOT. ok_veget) THEN
447     DO j = 1, knon
448     i = ni(j)
449     yqsol(j) = qsol(i)
450 guez 38 END DO
451 guez 62 ELSE
452     yqsol = 0.
453     END IF
454 guez 3
455 guez 62 DO k = 1, nsoilmx
456     DO j = 1, knon
457     i = ni(j)
458     ytsoil(j, k) = ftsoil(i, k, nsrf)
459 guez 38 END DO
460 guez 47 END DO
461 guez 3
462 guez 38 DO k = 1, klev
463     DO j = 1, knon
464     i = ni(j)
465 guez 62 ypaprs(j, k) = paprs(i, k)
466     ypplay(j, k) = pplay(i, k)
467     ydelp(j, k) = delp(i, k)
468     yu(j, k) = u(i, k)
469     yv(j, k) = v(i, k)
470     yt(j, k) = t(i, k)
471     yq(j, k) = q(i, k)
472 guez 38 END DO
473     END DO
474 guez 3
475 guez 62 ! calculer Cdrag et les coefficients d'echange
476     CALL coefkz(nsrf, knon, ypaprs, ypplay, ksta, ksta_ter, yts, yrugos, &
477     yu, yv, yt, yq, yqsurf, coefm(:knon, :), coefh(:knon, :))
478     IF (iflag_pbl == 1) THEN
479     CALL coefkz2(nsrf, knon, ypaprs, ypplay, yt, ycoefm0, ycoefh0)
480     coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
481     coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
482     END IF
483 guez 3
484 guez 70 ! on met un seuil pour coefm et coefh
485 guez 62 IF (nsrf == is_oce) THEN
486     coefm(:knon, 1) = min(coefm(:knon, 1), cdmmax)
487     coefh(:knon, 1) = min(coefh(:knon, 1), cdhmax)
488 guez 38 END IF
489 guez 3
490 guez 62 IF (ok_kzmin) THEN
491     ! Calcul d'une diffusion minimale pour les conditions tres stables
492     CALL coefkzmin(knon, ypaprs, ypplay, yu, yv, yt, yq, &
493 guez 70 coefm(:knon, 1), ycoefm0, ycoefh0)
494 guez 62 coefm(:knon, :) = max(coefm(:knon, :), ycoefm0(:knon, :))
495     coefh(:knon, :) = max(coefh(:knon, :), ycoefh0(:knon, :))
496 guez 98 END IF
497 guez 3
498 guez 62 IF (iflag_pbl >= 3) THEN
499     ! MELLOR ET YAMADA adapté à Mars, Richard Fournier et
500     ! Frédéric Hourdin
501     yzlay(:knon, 1) = rd * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
502     + ypplay(:knon, 1))) &
503     * (ypaprs(:knon, 1) - ypplay(:knon, 1)) / rg
504     DO k = 2, klev
505     yzlay(1:knon, k) = yzlay(1:knon, k-1) &
506     + rd * 0.5 * (yt(1:knon, k-1) + yt(1:knon, k)) &
507     / ypaprs(1:knon, k) &
508     * (ypplay(1:knon, k-1) - ypplay(1:knon, k)) / rg
509     END DO
510     DO k = 1, klev
511     yteta(1:knon, k) = yt(1:knon, k)*(ypaprs(1:knon, 1) &
512     / ypplay(1:knon, k))**rkappa * (1.+0.61*yq(1:knon, k))
513     END DO
514     yzlev(1:knon, 1) = 0.
515     yzlev(:knon, klev+1) = 2. * yzlay(:knon, klev) &
516     - yzlay(:knon, klev - 1)
517     DO k = 2, klev
518     yzlev(1:knon, k) = 0.5*(yzlay(1:knon, k)+yzlay(1:knon, k-1))
519     END DO
520     DO k = 1, klev + 1
521     DO j = 1, knon
522     i = ni(j)
523     yq2(j, k) = q2(i, k, nsrf)
524     END DO
525     END DO
526    
527     CALL ustarhb(knon, yu, yv, coefm(:knon, 1), yustar)
528    
529     IF (prt_level > 9) THEN
530     PRINT *, 'USTAR = ', yustar
531     END IF
532    
533     ! iflag_pbl peut être utilisé comme longueur de mélange
534    
535     IF (iflag_pbl >= 11) THEN
536     CALL vdif_kcay(knon, dtime, rg, rd, ypaprs, yt, yzlev, yzlay, &
537     yu, yv, yteta, coefm(:knon, 1), yq2, q2diag, ykmm, ykmn, &
538     yustar, iflag_pbl)
539     ELSE
540     CALL yamada4(knon, dtime, rg, yzlev, yzlay, yu, yv, yteta, &
541     coefm(:knon, 1), yq2, ykmm, ykmn, ykmq, yustar, iflag_pbl)
542     END IF
543    
544     coefm(:knon, 2:) = ykmm(:knon, 2:klev)
545     coefh(:knon, 2:) = ykmn(:knon, 2:klev)
546 guez 38 END IF
547 guez 3
548 guez 62 ! calculer la diffusion des vitesses "u" et "v"
549 guez 70 CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yu, ypaprs, &
550     ypplay, ydelp, y_d_u, y_flux_u)
551     CALL clvent(knon, dtime, yu1, yv1, coefm(:knon, :), yt, yv, ypaprs, &
552     ypplay, ydelp, y_d_v, y_flux_v)
553 guez 3
554 guez 62 ! pour le couplage
555     ytaux = y_flux_u(:, 1)
556     ytauy = y_flux_v(:, 1)
557 guez 3
558 guez 62 ! calculer la diffusion de "q" et de "h"
559 guez 72 CALL clqh(dtime, itap, jour, debut, rlat, knon, nsrf, ni, pctsrf, &
560     soil_model, ytsoil, yqsol, ok_veget, ocean, rmu0, co2_ppm, &
561     yrugos, yrugoro, yu1, yv1, coefh(:knon, :), yt, yq, yts, &
562     ypaprs, ypplay, ydelp, yrads, yalb, yalblw, ysnow, yqsurf, &
563     yrain_f, ysnow_f, yfder, ysolsw, yfluxlat, pctsrf_new, &
564     yagesno, y_d_t, y_d_q, y_d_ts, yz0_new, y_flux_t, y_flux_q, &
565     y_dflux_t, y_dflux_q, y_fqcalving, y_ffonte, y_run_off_lic_0, &
566     y_flux_o, y_flux_g, ytslab, y_seaice)
567 guez 3
568 guez 62 ! calculer la longueur de rugosite sur ocean
569     yrugm = 0.
570     IF (nsrf == is_oce) THEN
571     DO j = 1, knon
572     yrugm(j) = 0.018*coefm(j, 1)*(yu1(j)**2+yv1(j)**2)/rg + &
573     0.11*14E-6/sqrt(coefm(j, 1)*(yu1(j)**2+yv1(j)**2))
574     yrugm(j) = max(1.5E-05, yrugm(j))
575     END DO
576     END IF
577 guez 38 DO j = 1, knon
578 guez 62 y_dflux_t(j) = y_dflux_t(j)*ypct(j)
579     y_dflux_q(j) = y_dflux_q(j)*ypct(j)
580     yu1(j) = yu1(j)*ypct(j)
581     yv1(j) = yv1(j)*ypct(j)
582 guez 38 END DO
583 guez 3
584 guez 62 DO k = 1, klev
585     DO j = 1, knon
586     i = ni(j)
587     coefh(j, k) = coefh(j, k)*ypct(j)
588     coefm(j, k) = coefm(j, k)*ypct(j)
589     y_d_t(j, k) = y_d_t(j, k)*ypct(j)
590     y_d_q(j, k) = y_d_q(j, k)*ypct(j)
591     flux_t(i, k, nsrf) = y_flux_t(j, k)
592     flux_q(i, k, nsrf) = y_flux_q(j, k)
593     flux_u(i, k, nsrf) = y_flux_u(j, k)
594     flux_v(i, k, nsrf) = y_flux_v(j, k)
595     y_d_u(j, k) = y_d_u(j, k)*ypct(j)
596     y_d_v(j, k) = y_d_v(j, k)*ypct(j)
597     END DO
598 guez 38 END DO
599 guez 3
600 guez 62 evap(:, nsrf) = -flux_q(:, 1, nsrf)
601 guez 15
602 guez 62 albe(:, nsrf) = 0.
603     alblw(:, nsrf) = 0.
604     snow(:, nsrf) = 0.
605     qsurf(:, nsrf) = 0.
606     rugos(:, nsrf) = 0.
607     fluxlat(:, nsrf) = 0.
608 guez 38 DO j = 1, knon
609     i = ni(j)
610 guez 62 d_ts(i, nsrf) = y_d_ts(j)
611     albe(i, nsrf) = yalb(j)
612     alblw(i, nsrf) = yalblw(j)
613     snow(i, nsrf) = ysnow(j)
614     qsurf(i, nsrf) = yqsurf(j)
615     rugos(i, nsrf) = yz0_new(j)
616     fluxlat(i, nsrf) = yfluxlat(j)
617     IF (nsrf == is_oce) THEN
618     rugmer(i) = yrugm(j)
619     rugos(i, nsrf) = yrugm(j)
620     END IF
621     agesno(i, nsrf) = yagesno(j)
622     fqcalving(i, nsrf) = y_fqcalving(j)
623     ffonte(i, nsrf) = y_ffonte(j)
624     cdragh(i) = cdragh(i) + coefh(j, 1)
625     cdragm(i) = cdragm(i) + coefm(j, 1)
626     dflux_t(i) = dflux_t(i) + y_dflux_t(j)
627     dflux_q(i) = dflux_q(i) + y_dflux_q(j)
628     zu1(i) = zu1(i) + yu1(j)
629     zv1(i) = zv1(i) + yv1(j)
630 guez 38 END DO
631 guez 62 IF (nsrf == is_ter) THEN
632     DO j = 1, knon
633     i = ni(j)
634     qsol(i) = yqsol(j)
635     END DO
636     END IF
637     IF (nsrf == is_lic) THEN
638     DO j = 1, knon
639     i = ni(j)
640     run_off_lic_0(i) = y_run_off_lic_0(j)
641     END DO
642     END IF
643     !$$$ PB ajout pour soil
644     ftsoil(:, :, nsrf) = 0.
645     DO k = 1, nsoilmx
646     DO j = 1, knon
647     i = ni(j)
648     ftsoil(i, k, nsrf) = ytsoil(j, k)
649     END DO
650     END DO
651    
652 guez 38 DO j = 1, knon
653     i = ni(j)
654 guez 62 DO k = 1, klev
655     d_t(i, k) = d_t(i, k) + y_d_t(j, k)
656     d_q(i, k) = d_q(i, k) + y_d_q(j, k)
657     d_u(i, k) = d_u(i, k) + y_d_u(j, k)
658     d_v(i, k) = d_v(i, k) + y_d_v(j, k)
659 guez 70 ycoefh(i, k) = ycoefh(i, k) + coefh(j, k)
660 guez 62 END DO
661 guez 38 END DO
662 guez 62
663     !cc diagnostic t, q a 2m et u, v a 10m
664    
665 guez 38 DO j = 1, knon
666     i = ni(j)
667 guez 62 uzon(j) = yu(j, 1) + y_d_u(j, 1)
668     vmer(j) = yv(j, 1) + y_d_v(j, 1)
669     tair1(j) = yt(j, 1) + y_d_t(j, 1)
670     qair1(j) = yq(j, 1) + y_d_q(j, 1)
671     zgeo1(j) = rd*tair1(j)/(0.5*(ypaprs(j, 1)+ypplay(j, &
672     1)))*(ypaprs(j, 1)-ypplay(j, 1))
673     tairsol(j) = yts(j) + y_d_ts(j)
674     rugo1(j) = yrugos(j)
675     IF (nsrf == is_oce) THEN
676     rugo1(j) = rugos(i, nsrf)
677     END IF
678     psfce(j) = ypaprs(j, 1)
679     patm(j) = ypplay(j, 1)
680 guez 15
681 guez 62 qairsol(j) = yqsurf(j)
682 guez 38 END DO
683 guez 15
684 guez 62 CALL stdlevvar(klon, knon, nsrf, zxli, uzon, vmer, tair1, qair1, &
685     zgeo1, tairsol, qairsol, rugo1, psfce, patm, yt2m, yq2m, &
686     yt10m, yq10m, yu10m, yustar)
687 guez 3
688 guez 62 DO j = 1, knon
689     i = ni(j)
690     t2m(i, nsrf) = yt2m(j)
691     q2m(i, nsrf) = yq2m(j)
692 guez 3
693 guez 62 ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
694     u10m(i, nsrf) = (yu10m(j)*uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
695     v10m(i, nsrf) = (yu10m(j)*vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
696 guez 15
697 guez 62 END DO
698 guez 15
699 guez 62 CALL hbtm(knon, ypaprs, ypplay, yt2m, yt10m, yq2m, yq10m, yustar, &
700     y_flux_t, y_flux_q, yu, yv, yt, yq, ypblh, ycapcl, yoliqcl, &
701     ycteicl, ypblt, ytherm, ytrmb1, ytrmb2, ytrmb3, ylcl)
702 guez 15
703 guez 38 DO j = 1, knon
704     i = ni(j)
705 guez 62 pblh(i, nsrf) = ypblh(j)
706     plcl(i, nsrf) = ylcl(j)
707     capcl(i, nsrf) = ycapcl(j)
708     oliqcl(i, nsrf) = yoliqcl(j)
709     cteicl(i, nsrf) = ycteicl(j)
710     pblt(i, nsrf) = ypblt(j)
711     therm(i, nsrf) = ytherm(j)
712     trmb1(i, nsrf) = ytrmb1(j)
713     trmb2(i, nsrf) = ytrmb2(j)
714     trmb3(i, nsrf) = ytrmb3(j)
715 guez 38 END DO
716 guez 3
717 guez 38 DO j = 1, knon
718 guez 62 DO k = 1, klev + 1
719     i = ni(j)
720     q2(i, k, nsrf) = yq2(j, k)
721     END DO
722 guez 38 END DO
723 guez 62 !IM "slab" ocean
724 guez 47 IF (nsrf == is_oce) THEN
725 guez 62 DO j = 1, knon
726     ! on projette sur la grille globale
727     i = ni(j)
728     IF (pctsrf_new(i, is_oce)>epsfra) THEN
729     flux_o(i) = y_flux_o(j)
730     ELSE
731     flux_o(i) = 0.
732     END IF
733     END DO
734 guez 38 END IF
735 guez 62
736     IF (nsrf == is_sic) THEN
737     DO j = 1, knon
738     i = ni(j)
739     ! On pondère lorsque l'on fait le bilan au sol :
740     IF (pctsrf_new(i, is_sic)>epsfra) THEN
741     flux_g(i) = y_flux_g(j)
742     ELSE
743     flux_g(i) = 0.
744     END IF
745     END DO
746    
747     END IF
748     IF (ocean == 'slab ') THEN
749     IF (nsrf == is_oce) THEN
750     tslab(1:klon) = ytslab(1:klon)
751     seaice(1:klon) = y_seaice(1:klon)
752     END IF
753     END IF
754     end IF if_knon
755 guez 49 END DO loop_surface
756 guez 15
757 guez 38 ! On utilise les nouvelles surfaces
758 guez 15
759 guez 38 rugos(:, is_oce) = rugmer
760     pctsrf = pctsrf_new
761 guez 15
762 guez 38 END SUBROUTINE clmain
763 guez 15
764 guez 38 end module clmain_m

  ViewVC Help
Powered by ViewVC 1.1.21