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

Annotation of /trunk/phylmd/clmain.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 25 - (hide annotations)
Fri Mar 5 16:43:45 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/phylmd/clmain.f90
File size: 27774 byte(s)
Simplified "etat0_lim.sh" and "gcm.sh" because the full versions
depended on personal arrangements for directories and machines.

Translated included files into modules. Encapsulated procedures into modules.

Moved variables from module "comgeom" to local variables of
"inigeom". Deleted some unused variables in "comgeom".

Moved variable "day_ini" from module "temps" to module "dynetat0_m".

Removed useless test on variable "time" and useless "close" statement
in procedure "leapfrog".

Removed useless call to "inigeom" in procedure "limit".

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

  ViewVC Help
Powered by ViewVC 1.1.21