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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (show annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years ago) by guez
Original Path: trunk/libf/phylmd/clmain.f
File size: 29420 byte(s)
Added NetCDF directory "/home/guez/include" in "g95.mk" and
"nag_tools.mk".

Added some "intent" attributes in "PVtheta", "advtrac", "caladvtrac",
"calfis", "diagedyn", "dissip", "vlspltqs", "aeropt", "ajsec",
"calltherm", "clmain", "cltrac", "cltracrn", "concvl", "conema3",
"conflx", "fisrtilp", "newmicro", "nuage", "diagcld1", "diagcld2",
"drag_noro", "lift_noro", "SUGWD", "physiq", "phytrac", "radlwsw", "thermcell".

Removed the case "ierr == 0" in "abort_gcm"; moved call to "histclo"
and messages for end of run from "abort_gcm" to "gcm"; replaced call
to "abort_gcm" in "leapfrog" by exit from outer loop.

In "calfis": removed argument "pp" and variable "unskap"; changed
"pksurcp" from scalar to rank 2; use "pressure_var"; rewrote
computation of "zplev", "zplay", "ztfi", "pcvgt" using "dyn_phy";
added computation of "pls".

Removed unused variable in "dynredem0".

In "exner_hyb": changed "dellta" from scalar to rank 1; replaced call
to "ssum" by call to "sum"; removed variables "xpn" and "xps";
replaced some loops by array expressions.

In "leapfrog": use "pressure_var"; deleted variables "p", "longcles".

Converted common blocks "YOECUMF", "YOEGWD" to modules.

Removed argument "pplay" in "cvltr", "diagetpq", "nflxtr".

Created module "raddimlw" from include file "raddimlw.h".

Corrected call to "new_unit" in "test_disvert".

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 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),rugoro(klon)
145 REAL cdragh(klon), cdragm(klon)
146 integer jour ! jour de l'annee en cours
147 real rmu0(klon) ! cosinus de l'angle solaire zenithal
148 REAL co2_ppm ! taux CO2 atmosphere
149 LOGICAL, intent(in):: debut
150 logical, intent(in):: lafin
151 logical ok_veget
152 character*6 ocean
153 integer npas, nexca
154 c
155 REAL pctsrf(klon,nbsrf)
156 REAL ts(klon,nbsrf)
157 REAL d_ts(klon,nbsrf)
158 REAL snow(klon,nbsrf)
159 REAL qsurf(klon,nbsrf)
160 REAL evap(klon,nbsrf)
161 REAL albe(klon,nbsrf)
162 REAL alblw(klon,nbsrf)
163 c$$$ PB
164 REAL fluxlat(klon,nbsrf)
165 C
166 real rain_f(klon), snow_f(klon)
167 REAL fder(klon)
168 cIM cf. JLD REAL sollw(klon), solsw(klon), sollwdown(klon)
169 REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)
170 REAL rugos(klon,nbsrf)
171 C la nouvelle repartition des surfaces sortie de l'interface
172 REAL pctsrf_new(klon,nbsrf)
173 cAA
174 REAL zcoefh(klon,klev)
175 REAL zu1(klon)
176 REAL zv1(klon)
177 cAA
178 c$$$ PB ajout pour soil
179 LOGICAL soil_model
180 cIM ajout seuils cdrm, cdrh
181 REAL cdmmax, cdhmax
182 cIM: 261103
183 REAL ksta, ksta_ter
184 LOGICAL ok_kzmin
185 cIM: 261103
186 REAL ftsoil(klon,nsoilmx,nbsrf)
187 REAL ytsoil(klon,nsoilmx)
188 REAL qsol(klon)
189 c======================================================================
190 EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
191 c======================================================================
192 REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
193 REAL yalb(klon)
194 REAL yalblw(klon)
195 REAL yu1(klon), yv1(klon)
196 real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
197 real yrain_f(klon), ysnow_f(klon)
198 real ysollw(klon), ysolsw(klon), ysollwdown(klon)
199 real yfder(klon), ytaux(klon), ytauy(klon)
200 REAL yrugm(klon), yrads(klon),yrugoro(klon)
201 c$$$ PB
202 REAL yfluxlat(klon)
203 C
204 REAL y_d_ts(klon)
205 REAL y_d_t(klon, klev), y_d_q(klon, klev)
206 REAL y_d_u(klon, klev), y_d_v(klon, klev)
207 REAL y_flux_t(klon,klev), y_flux_q(klon,klev)
208 REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
209 REAL y_dflux_t(klon), y_dflux_q(klon)
210 REAL ycoefh(klon,klev), ycoefm(klon,klev)
211 REAL yu(klon,klev), yv(klon,klev)
212 REAL yt(klon,klev), yq(klon,klev)
213 REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
214 c
215 LOGICAL ok_nonloc
216 PARAMETER (ok_nonloc=.FALSE.)
217 REAL ycoefm0(klon,klev), ycoefh0(klon,klev)
218
219 cIM 081204 hcl_Anne ? BEG
220 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 real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)
224 real q2diag(klon,klev+1)
225 cIM 081204 real yustar(klon),y_cd_m(klon),y_cd_h(klon)
226 cIM 081204 hcl_Anne ? END
227 c
228 REAL u1lay(klon), v1lay(klon)
229 REAL delp(klon,klev)
230 INTEGER i, k, nsrf
231 cAA INTEGER it
232 INTEGER ni(klon), knon, j
233 c Introduction d'une variable "pourcentage potentiel" pour tenir compte
234 c des eventuelles apparitions et/ou disparitions de la glace de mer
235 REAL pctsrf_pot(klon,nbsrf)
236
237 c======================================================================
238 REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
239 c======================================================================
240 c
241 c maf pour sorties IOISPL en cas de debugagage
242 c
243 CHARACTER*80 cldebug
244 SAVE cldebug
245 CHARACTER*8 cl_surf(nbsrf)
246 SAVE cl_surf
247 INTEGER nhoridbg, nidbg
248 SAVE nhoridbg, nidbg
249 INTEGER ndexbg(iim*(jjm+1))
250 REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
251 REAL tabindx(klon)
252 REAL debugtab(iim,jjm+1)
253 LOGICAL first_appel
254 SAVE first_appel
255 DATA first_appel/.true./
256 LOGICAL debugindex
257 SAVE debugindex
258 DATA debugindex/.false./
259 integer idayref
260 REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
261 REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
262 c
263 REAL yt2m(klon), yq2m(klon), yu10m(klon)
264 REAL yustar(klon)
265 c -- LOOP
266 REAL yu10mx(klon)
267 REAL yu10my(klon)
268 REAL ywindsp(klon)
269 c -- LOOP
270 c
271 REAL yt10m(klon), yq10m(klon)
272 cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui
273 c permet de sortir les grdeurs par sous surface)
274 REAL pblh(klon,nbsrf)
275 REAL plcl(klon,nbsrf)
276 REAL capCL(klon,nbsrf)
277 REAL oliqCL(klon,nbsrf)
278 REAL cteiCL(klon,nbsrf)
279 REAL pblT(klon,nbsrf)
280 REAL therm(klon,nbsrf)
281 REAL trmb1(klon,nbsrf)
282 REAL trmb2(klon,nbsrf)
283 REAL trmb3(klon,nbsrf)
284 REAL ypblh(klon)
285 REAL ylcl(klon)
286 REAL ycapCL(klon)
287 REAL yoliqCL(klon)
288 REAL ycteiCL(klon)
289 REAL ypblT(klon)
290 REAL ytherm(klon)
291 REAL ytrmb1(klon)
292 REAL ytrmb2(klon)
293 REAL ytrmb3(klon)
294 REAL y_cd_h(klon), y_cd_m(klon)
295 c REAL ygamt(klon,2:klev) ! contre-gradient pour temperature
296 c REAL ygamq(klon,2:klev) ! contre-gradient pour humidite
297 REAL uzon(klon), vmer(klon)
298 REAL tair1(klon), qair1(klon), tairsol(klon)
299 REAL psfce(klon), patm(klon)
300 c
301 REAL qairsol(klon), zgeo1(klon)
302 REAL rugo1(klon)
303 c
304 LOGICAL zxli ! utiliser un jeu de fonctions simples
305 PARAMETER (zxli=.FALSE.)
306 c
307 REAL zt, zqs, zdelta, zcor
308 REAL t_coup
309 PARAMETER(t_coup=273.15)
310 C
311 character (len = 20) :: modname = 'clmain'
312 LOGICAL check
313 PARAMETER (check=.false.)
314
315
316 c initialisation Anne
317 ytherm(:) = 0.
318 C
319 if (check) THEN
320 write(*,*) modname,' klon=',klon
321 CC call flush(6)
322 endif
323 IF (debugindex .and. first_appel) THEN
324 first_appel=.false.
325 !
326 ! initialisation sorties netcdf
327 !
328 idayref = day_ini
329 CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
330 CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
331 DO i = 1, iim
332 zx_lon(i,1) = rlon(i+1)
333 zx_lon(i,jjm+1) = rlon(i+1)
334 ENDDO
335 CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
336 cldebug='sous_index'
337 CALL histbeg_totreg(cldebug, iim,zx_lon(:,1),jjm+1,
338 $ zx_lat(1,:),1,iim,1,jjm
339 $ +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)
340 ! no vertical axis
341 cl_surf(1)='ter'
342 cl_surf(2)='lic'
343 cl_surf(3)='oce'
344 cl_surf(4)='sic'
345 DO nsrf=1,nbsrf
346 CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,
347 $ jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
348 END DO
349 CALL histend(nidbg)
350 CALL histsync(nidbg)
351 ENDIF
352
353 DO k = 1, klev ! epaisseur de couche
354 DO i = 1, klon
355 delp(i,k) = paprs(i,k)-paprs(i,k+1)
356 ENDDO
357 ENDDO
358 DO i = 1, klon ! vent de la premiere couche
359 zx_alf1 = 1.0
360 zx_alf2 = 1.0 - zx_alf1
361 u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
362 v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
363 ENDDO
364 c
365 c initialisation:
366 c
367 DO i = 1, klon
368 rugmer(i) = 0.0
369 cdragh(i) = 0.0
370 cdragm(i) = 0.0
371 dflux_t(i) = 0.0
372 dflux_q(i) = 0.0
373 zu1(i) = 0.0
374 zv1(i) = 0.0
375 ENDDO
376 ypct = 0.0
377 yts = 0.0
378 ysnow = 0.0
379 yqsurf = 0.0
380 yalb = 0.0
381 yalblw = 0.0
382 yrain_f = 0.0
383 ysnow_f = 0.0
384 yfder = 0.0
385 ytaux = 0.0
386 ytauy = 0.0
387 ysolsw = 0.0
388 ysollw = 0.0
389 ysollwdown = 0.0
390 yrugos = 0.0
391 yu1 = 0.0
392 yv1 = 0.0
393 yrads = 0.0
394 ypaprs = 0.0
395 ypplay = 0.0
396 ydelp = 0.0
397 yu = 0.0
398 yv = 0.0
399 yt = 0.0
400 yq = 0.0
401 pctsrf_new = 0.0
402 y_flux_u = 0.0
403 y_flux_v = 0.0
404 C$$ PB
405 y_dflux_t = 0.0
406 y_dflux_q = 0.0
407 ytsoil = 999999.
408 yrugoro = 0.
409 c -- LOOP
410 yu10mx = 0.0
411 yu10my = 0.0
412 ywindsp = 0.0
413 c -- LOOP
414 DO nsrf = 1, nbsrf
415 DO i = 1, klon
416 d_ts(i,nsrf) = 0.0
417 ENDDO
418 END DO
419 C§§§ PB
420 yfluxlat=0.
421 flux_t = 0.
422 flux_q = 0.
423 flux_u = 0.
424 flux_v = 0.
425 DO k = 1, klev
426 DO i = 1, klon
427 d_t(i,k) = 0.0
428 d_q(i,k) = 0.0
429 c$$$ flux_t(i,k) = 0.0
430 c$$$ flux_q(i,k) = 0.0
431 d_u(i,k) = 0.0
432 d_v(i,k) = 0.0
433 c$$$ flux_u(i,k) = 0.0
434 c$$$ flux_v(i,k) = 0.0
435 zcoefh(i,k) = 0.0
436 ENDDO
437 ENDDO
438 cAA IF (itr.GE.1) THEN
439 cAA DO it = 1, itr
440 cAA DO k = 1, klev
441 cAA DO i = 1, klon
442 cAA d_tr(i,k,it) = 0.0
443 cAA ENDDO
444 cAA ENDDO
445 cAA ENDDO
446 cAA ENDIF
447
448 c
449 c Boucler sur toutes les sous-fractions du sol:
450 c
451 C Initialisation des "pourcentages potentiels". On considere ici qu'on
452 C peut avoir potentiellementdela glace sur tout le domaine oceanique
453 C (a affiner)
454
455 pctsrf_pot = pctsrf
456 pctsrf_pot(:,is_oce) = 1. - zmasq(:)
457 pctsrf_pot(:,is_sic) = 1. - zmasq(:)
458
459 DO 99999 nsrf = 1, nbsrf
460
461 c chercher les indices:
462 DO j = 1, klon
463 ni(j) = 0
464 ENDDO
465 knon = 0
466 DO i = 1, klon
467
468 C pour determiner le domaine a traiter on utilise les surfaces "potentielles"
469 C
470 IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
471 knon = knon + 1
472 ni(knon) = i
473 ENDIF
474 ENDDO
475 c
476 if (check) THEN
477 write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon
478 CC call flush(6)
479 endif
480 c
481 c variables pour avoir une sortie IOIPSL des INDEX
482 c
483 IF (debugindex) THEN
484 tabindx(:)=0.
485 c tabindx(1:knon)=(/FLOAT(i),i=1:knon/)
486 DO i=1,knon
487 tabindx(1:knon)=FLOAT(i)
488 END DO
489 debugtab(:,:)=0.
490 ndexbg(:)=0
491 CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)
492 CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)
493 $ ,ndexbg)
494 ENDIF
495 IF (knon.EQ.0) GOTO 99999
496 DO j = 1, knon
497 i = ni(j)
498 ypct(j) = pctsrf(i,nsrf)
499 yts(j) = ts(i,nsrf)
500 cIM "slab" ocean
501 c PRINT *, 'tslab = ', i, tslab(i)
502 ytslab(i) = tslab(i)
503 c
504 ysnow(j) = snow(i,nsrf)
505 yqsurf(j) = qsurf(i,nsrf)
506 yalb(j) = albe(i,nsrf)
507 yalblw(j) = alblw(i,nsrf)
508 yrain_f(j) = rain_f(i)
509 ysnow_f(j) = snow_f(i)
510 yagesno(j) = agesno(i,nsrf)
511 yfder(j) = fder(i)
512 ytaux(j) = flux_u(i,1,nsrf)
513 ytauy(j) = flux_v(i,1,nsrf)
514 ysolsw(j) = solsw(i,nsrf)
515 ysollw(j) = sollw(i,nsrf)
516 ysollwdown(j) = sollwdown(i)
517 yrugos(j) = rugos(i,nsrf)
518 yrugoro(j) = rugoro(i)
519 yu1(j) = u1lay(i)
520 yv1(j) = v1lay(i)
521 yrads(j) = ysolsw(j)+ ysollw(j)
522 ypaprs(j,klev+1) = paprs(i,klev+1)
523 y_run_off_lic_0(j) = run_off_lic_0(i)
524 c -- LOOP
525 yu10mx(j) = u10m(i,nsrf)
526 yu10my(j) = v10m(i,nsrf)
527 ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )
528 c -- LOOP
529 END DO
530 C
531 C IF bucket model for continent, copy soil water content
532 IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN
533 DO j = 1, knon
534 i = ni(j)
535 yqsol(j) = qsol(i)
536 END DO
537 ELSE
538 yqsol(:)=0.
539 ENDIF
540 c$$$ PB ajour pour soil
541 DO k = 1, nsoilmx
542 DO j = 1, knon
543 i = ni(j)
544 ytsoil(j,k) = ftsoil(i,k,nsrf)
545 END DO
546 END DO
547 DO k = 1, klev
548 DO j = 1, knon
549 i = ni(j)
550 ypaprs(j,k) = paprs(i,k)
551 ypplay(j,k) = pplay(i,k)
552 ydelp(j,k) = delp(i,k)
553 yu(j,k) = u(i,k)
554 yv(j,k) = v(i,k)
555 yt(j,k) = t(i,k)
556 yq(j,k) = q(i,k)
557 ENDDO
558 ENDDO
559 c
560 c
561 c calculer Cdrag et les coefficients d'echange
562 CALL coefkz(nsrf, knon, ypaprs, ypplay,
563 cIM 261103
564 . ksta, ksta_ter,
565 cIM 261103
566 . yts, yrugos, yu, yv, yt, yq,
567 . yqsurf,
568 . ycoefm, ycoefh)
569 cIM 081204 BEG
570 cCR test
571 if (iflag_pbl.eq.1) then
572 cIM 081204 END
573 CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,
574 . ycoefm0, ycoefh0)
575 DO k = 1, klev
576 DO i = 1, knon
577 ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
578 ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
579 ENDDO
580 ENDDO
581 endif
582 c
583 cIM cf JLD : on seuille ycoefm et ycoefh
584 if (nsrf.eq.is_oce) then
585 do j=1,knon
586 c ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)
587 ycoefm(j,1)=min(ycoefm(j,1),cdmmax)
588 c ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)
589 ycoefh(j,1)=min(ycoefh(j,1),cdhmax)
590 enddo
591 endif
592
593 c
594 cIM: 261103
595 if (ok_kzmin) THEN
596 cIM cf FH: 201103 BEG
597 c Calcul d'une diffusion minimale pour les conditions tres stables.
598 call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm
599 . ,ycoefm0,ycoefh0)
600 c call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ ')
601 c call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN ')
602
603 if ( 1.eq.1 ) then
604 DO k = 1, klev
605 DO i = 1, knon
606 ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
607 ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
608 ENDDO
609 ENDDO
610 endif
611 cIM cf FH: 201103 END
612 endif !ok_kzmin
613 cIM: 261103
614
615
616 IF (iflag_pbl.ge.3) then
617
618 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
619 c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
620 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
621
622 yzlay(1:knon,1)=
623 . RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1)))
624 . *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
625 do k=2,klev
626 yzlay(1:knon,k)=
627 . yzlay(1:knon,k-1)+RD*0.5*(yt(1:knon,k-1)+yt(1:knon,k))
628 . /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG
629 enddo
630 do k=1,klev
631 yteta(1:knon,k)=
632 . yt(1:knon,k)*(ypaprs(1:knon,1)/ypplay(1:knon,k))**rkappa
633 . *(1.+0.61*yq(1:knon,k))
634 enddo
635 yzlev(1:knon,1)=0.
636 yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
637 do k=2,klev
638 yzlev(1:knon,k)=0.5*(yzlay(1:knon,k)+yzlay(1:knon,k-1))
639 enddo
640 DO k = 1, klev+1
641 DO j = 1, knon
642 i = ni(j)
643 yq2(j,k)=q2(i,k,nsrf)
644 enddo
645 enddo
646
647
648 c Bug introduit volontairement pour converger avec les resultats
649 c du papier sur les thermiques.
650 if (1.eq.1) then
651 y_cd_m(1:knon) = ycoefm(1:knon,1)
652 y_cd_h(1:knon) = ycoefh(1:knon,1)
653 else
654 y_cd_h(1:knon) = ycoefm(1:knon,1)
655 y_cd_m(1:knon) = ycoefh(1:knon,1)
656 endif
657 call ustarhb(knon,yu,yv,y_cd_m, yustar)
658
659 if (prt_level > 9) THEN
660 WRITE(lunout,*)'USTAR = ',yustar
661 ENDIF
662
663 c iflag_pbl peut etre utilise comme longuer de melange
664
665 if (iflag_pbl.ge.11) then
666 call vdif_kcay(knon,dtime,rg,rd,ypaprs,yt
667 s ,yzlev,yzlay,yu,yv,yteta
668 s ,y_cd_m,yq2,q2diag,ykmm,ykmn,yustar,
669 s iflag_pbl)
670 else
671 call yamada4(knon,dtime,rg,rd,ypaprs,yt
672 s ,yzlev,yzlay,yu,yv,yteta
673 s ,y_cd_m,yq2,ykmm,ykmn,ykmq,yustar,
674 s iflag_pbl)
675 endif
676
677 ycoefm(1:knon,1)=y_cd_m(1:knon)
678 ycoefh(1:knon,1)=y_cd_h(1:knon)
679 ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
680 ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
681
682
683 ENDIF
684
685 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
686 c calculer la diffusion des vitesses "u" et "v"
687 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
688
689 CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
690 s y_d_u,y_flux_u)
691 CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
692 s y_d_v,y_flux_v)
693
694 c pour le couplage
695 ytaux = y_flux_u(:,1)
696 ytauy = y_flux_v(:,1)
697
698 c FH modif sur le cdrag temperature
699 c$$$PB : déplace dans clcdrag
700 c$$$ do i=1,knon
701 c$$$ ycoefh(i,1)=ycoefm(i,1)*0.8
702 c$$$ enddo
703
704 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
705 c calculer la diffusion de "q" et de "h"
706 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
707 CALL clqh(dtime, itap, date0,jour, debut,lafin,
708 e rlon, rlat, cufi, cvfi,
709 e knon, nsrf, ni, pctsrf,
710 e soil_model, ytsoil,yqsol,
711 e ok_veget, ocean, npas, nexca,
712 e rmu0, co2_ppm, yrugos, yrugoro,
713 e yu1, yv1, ycoefh,
714 e yt,yq,yts,ypaprs,ypplay,
715 e ydelp,yrads,yalb, yalblw, ysnow, yqsurf,
716 e yrain_f, ysnow_f, yfder, ytaux, ytauy,
717 c -- LOOP
718 e ywindsp,
719 c -- LOOP
720 c$$$ e ysollw, ysolsw,
721 e ysollw, ysollwdown, ysolsw,yfluxlat,
722 s pctsrf_new, yagesno,
723 s y_d_t, y_d_q, y_d_ts, yz0_new,
724 s y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,
725 s y_fqcalving,y_ffonte,y_run_off_lic_0,
726 cIM "slab" ocean
727 s y_flux_o, y_flux_g, ytslab, y_seaice)
728 c
729 c calculer la longueur de rugosite sur ocean
730 yrugm=0.
731 IF (nsrf.EQ.is_oce) THEN
732 DO j = 1, knon
733 yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG
734 $ + 0.11*14e-6 / sqrt(ycoefm(j,1) * (yu1(j)**2+yv1(j)**2))
735 yrugm(j) = MAX(1.5e-05,yrugm(j))
736 ENDDO
737 ENDIF
738 DO j = 1, knon
739 y_dflux_t(j) = y_dflux_t(j) * ypct(j)
740 y_dflux_q(j) = y_dflux_q(j) * ypct(j)
741 yu1(j) = yu1(j) * ypct(j)
742 yv1(j) = yv1(j) * ypct(j)
743 ENDDO
744 c
745 DO k = 1, klev
746 DO j = 1, knon
747 i = ni(j)
748 ycoefh(j,k) = ycoefh(j,k) * ypct(j)
749 ycoefm(j,k) = ycoefm(j,k) * ypct(j)
750 y_d_t(j,k) = y_d_t(j,k) * ypct(j)
751 y_d_q(j,k) = y_d_q(j,k) * ypct(j)
752 C§§§ PB
753 flux_t(i,k,nsrf) = y_flux_t(j,k)
754 flux_q(i,k,nsrf) = y_flux_q(j,k)
755 flux_u(i,k,nsrf) = y_flux_u(j,k)
756 flux_v(i,k,nsrf) = y_flux_v(j,k)
757 c$$$ PB y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
758 c$$$ PB y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
759 y_d_u(j,k) = y_d_u(j,k) * ypct(j)
760 y_d_v(j,k) = y_d_v(j,k) * ypct(j)
761 c$$$ PB y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
762 c$$$ PB y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
763 ENDDO
764 ENDDO
765
766
767 evap(:,nsrf) = - flux_q(:,1,nsrf)
768 c
769 albe(:, nsrf) = 0.
770 alblw(:, nsrf) = 0.
771 snow(:, nsrf) = 0.
772 qsurf(:, nsrf) = 0.
773 rugos(:, nsrf) = 0.
774 fluxlat(:,nsrf) = 0.
775 DO j = 1, knon
776 i = ni(j)
777 d_ts(i,nsrf) = y_d_ts(j)
778 albe(i,nsrf) = yalb(j)
779 alblw(i,nsrf) = yalblw(j)
780 snow(i,nsrf) = ysnow(j)
781 qsurf(i,nsrf) = yqsurf(j)
782 rugos(i,nsrf) = yz0_new(j)
783 fluxlat(i,nsrf) = yfluxlat(j)
784 c$$$ pb rugmer(i) = yrugm(j)
785 IF (nsrf .EQ. is_oce) then
786 rugmer(i) = yrugm(j)
787 rugos(i,nsrf) = yrugm(j)
788 endif
789 cIM cf JLD ??
790 agesno(i,nsrf) = yagesno(j)
791 fqcalving(i,nsrf) = y_fqcalving(j)
792 ffonte(i,nsrf) = y_ffonte(j)
793 cdragh(i) = cdragh(i) + ycoefh(j,1)
794 cdragm(i) = cdragm(i) + ycoefm(j,1)
795 dflux_t(i) = dflux_t(i) + y_dflux_t(j)
796 dflux_q(i) = dflux_q(i) + y_dflux_q(j)
797 zu1(i) = zu1(i) + yu1(j)
798 zv1(i) = zv1(i) + yv1(j)
799 END DO
800 IF ( nsrf .eq. is_ter ) THEN
801 DO j = 1, knon
802 i = ni(j)
803 qsol(i) = yqsol(j)
804 END DO
805 END IF
806 IF ( nsrf .eq. is_lic ) THEN
807 DO j = 1, knon
808 i = ni(j)
809 run_off_lic_0(i) = y_run_off_lic_0(j)
810 END DO
811 END IF
812 c$$$ PB ajout pour soil
813 ftsoil(:,:,nsrf) = 0.
814 DO k = 1, nsoilmx
815 DO j = 1, knon
816 i = ni(j)
817 ftsoil(i, k, nsrf) = ytsoil(j,k)
818 END DO
819 END DO
820 c
821 DO j = 1, knon
822 i = ni(j)
823 DO k = 1, klev
824 d_t(i,k) = d_t(i,k) + y_d_t(j,k)
825 d_q(i,k) = d_q(i,k) + y_d_q(j,k)
826 c$$$ PB flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
827 c$$$ flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
828 d_u(i,k) = d_u(i,k) + y_d_u(j,k)
829 d_v(i,k) = d_v(i,k) + y_d_v(j,k)
830 c$$$ PB flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
831 c$$$ flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
832 zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
833 ENDDO
834 ENDDO
835 c
836 c
837 ccc diagnostic t,q a 2m et u, v a 10m
838 c
839 DO j=1, knon
840 i = ni(j)
841 uzon(j) = yu(j,1) + y_d_u(j,1)
842 vmer(j) = yv(j,1) + y_d_v(j,1)
843 tair1(j) = yt(j,1) + y_d_t(j,1)
844 qair1(j) = yq(j,1) + y_d_q(j,1)
845 zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1)))
846 & * (ypaprs(j,1)-ypplay(j,1))
847 tairsol(j) = yts(j) + y_d_ts(j)
848 rugo1(j) = yrugos(j)
849 IF(nsrf.EQ.is_oce) THEN
850 rugo1(j) = rugos(i,nsrf)
851 ENDIF
852 psfce(j)=ypaprs(j,1)
853 patm(j)=ypplay(j,1)
854 c
855 qairsol(j) = yqsurf(j)
856 ENDDO
857 c
858 CALL stdlevvar(klon, knon, nsrf, zxli,
859 & uzon, vmer, tair1, qair1, zgeo1,
860 & tairsol, qairsol, rugo1, psfce, patm,
861 cIM & yt2m, yq2m, yu10m)
862 & yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
863 cIM 081204 END
864 c
865 c
866 DO j=1, knon
867 i = ni(j)
868 t2m(i,nsrf)=yt2m(j)
869
870 c
871 q2m(i,nsrf)=yq2m(j)
872 c
873 c u10m, v10m : composantes du vent a 10m sans spirale de Ekman
874 u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
875 v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
876 c
877 ENDDO
878 c
879 cIM cf AM : pbl, HBTM
880 DO i = 1, knon
881 y_cd_h(i) = ycoefh(i,1)
882 y_cd_m(i) = ycoefm(i,1)
883 ENDDO
884 c print*,'appel hbtm2'
885 CALL HBTM(knon, ypaprs, ypplay,
886 . yt2m,yt10m,yq2m,yq10m,yustar,
887 . y_flux_t,y_flux_q,yu,yv,yt,yq,
888 . ypblh,ycapCL,yoliqCL,ycteiCL,ypblT,
889 . ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
890 c print*,'fin hbtm2'
891 c
892 DO j=1, knon
893 i = ni(j)
894 pblh(i,nsrf) = ypblh(j)
895 plcl(i,nsrf) = ylcl(j)
896 capCL(i,nsrf) = ycapCL(j)
897 oliqCL(i,nsrf) = yoliqCL(j)
898 cteiCL(i,nsrf) = ycteiCL(j)
899 pblT(i,nsrf) = ypblT(j)
900 therm(i,nsrf) = ytherm(j)
901 trmb1(i,nsrf) = ytrmb1(j)
902 trmb2(i,nsrf) = ytrmb2(j)
903 trmb3(i,nsrf) = ytrmb3(j)
904 ENDDO
905 c
906
907 do j=1,knon
908 do k=1,klev+1
909 i=ni(j)
910 q2(i,k,nsrf)=yq2(j,k)
911 enddo
912 enddo
913 cIM "slab" ocean
914 IF(OCEAN.EQ.'slab '.OR.OCEAN.EQ.'force ') THEN
915 IF (nsrf.EQ.is_oce) THEN
916 DO j = 1, knon
917 c on projette sur la grille globale
918 i = ni(j)
919 IF(pctsrf_new(i,is_oce).GT.epsfra) THEN
920 flux_o(i) = y_flux_o(j)
921 ELSE
922 flux_o(i) = 0.
923 ENDIF
924 ENDDO
925 ENDIF
926 c
927 IF (nsrf.EQ.is_sic) THEN
928 DO j = 1, knon
929 i = ni(j)
930 cIM 230604 on pondere lorsque l'on fait le bilan au sol : flux_g(i) = y_flux_g(j)*ypct(j)
931 IF(pctsrf_new(i,is_sic).GT.epsfra) THEN
932 flux_g(i) = y_flux_g(j)
933 ELSE
934 flux_g(i) = 0.
935 ENDIF
936 ENDDO
937 ENDIF !nsrf.EQ.is_sic
938 ENDIF !OCEAN
939 c
940 IF(OCEAN.EQ.'slab ') THEN
941 IF(nsrf.EQ.is_oce) then
942 tslab(1:klon) = ytslab(1:klon)
943 seaice(1:klon) = y_seaice(1:klon)
944 ENDIF !nsrf
945 ENDIF !OCEAN
946 99999 CONTINUE
947 C
948 C On utilise les nouvelles surfaces
949 C A rajouter: conservation de l'albedo
950 C
951 rugos(:,is_oce) = rugmer
952 pctsrf = pctsrf_new
953
954 RETURN
955 END

  ViewVC Help
Powered by ViewVC 1.1.21