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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21