/[lmdze]/trunk/libf/phylmd/clmain.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/clmain.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 14 - (hide 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 guez 3 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 guez 12 REAL, intent(in):: dtime
120 guez 3 real date0
121 guez 7 integer, intent(in):: itap
122 guez 3 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 guez 10 real, intent(in):: pplay(klon,klev)
127 guez 3 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 guez 13 REAL rugmer(klon), agesno(klon,nbsrf)
145     real, intent(in):: rugoro(klon)
146 guez 3 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 guez 12 character(len=*), intent(IN):: ocean
154 guez 3 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 guez 12 LOGICAL, intent(in):: soil_model
181 guez 3 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 guez 13 LOGICAL:: debugindex = .false.
258 guez 3 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 guez 14 DO nsrf = 1, nbsrf
459 guez 3
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 guez 14 IF (knon.EQ.0) cycle
495 guez 3 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 guez 12 print *,'USTAR = ',yustar
660 guez 3 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 guez 14 end do
944 guez 3 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