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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 70 - (hide annotations)
Mon Jun 24 15:39:52 2013 UTC (10 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/clmain.f90
File size: 25932 byte(s)
In procedure, "addfi" access directly the module variable "dtphys"
instead of going through an argument.

In "conflx", do not create a local variable for temperature with
reversed order of vertical levels. Instead, give an actual argument
with reversed order in "physiq".

Changed names of variables "rmd" and "rmv" from module "suphec_m" to
"md" and "mv".

In "hgardfou", print only the first temperature out of range found.

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

  ViewVC Help
Powered by ViewVC 1.1.21