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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (hide annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 2 months ago) by guez
Original Path: trunk/libf/phylmd/clmain.f
File size: 29462 byte(s)
Initial import
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     REAL dtime
120     real date0
121     integer 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 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     ccc zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
360     zx_alf1 = 1.0
361     zx_alf2 = 1.0 - zx_alf1
362     u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
363     v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
364     ENDDO
365     c
366     c initialisation:
367     c
368     DO i = 1, klon
369     rugmer(i) = 0.0
370     cdragh(i) = 0.0
371     cdragm(i) = 0.0
372     dflux_t(i) = 0.0
373     dflux_q(i) = 0.0
374     zu1(i) = 0.0
375     zv1(i) = 0.0
376     ENDDO
377     ypct = 0.0
378     yts = 0.0
379     ysnow = 0.0
380     yqsurf = 0.0
381     yalb = 0.0
382     yalblw = 0.0
383     yrain_f = 0.0
384     ysnow_f = 0.0
385     yfder = 0.0
386     ytaux = 0.0
387     ytauy = 0.0
388     ysolsw = 0.0
389     ysollw = 0.0
390     ysollwdown = 0.0
391     yrugos = 0.0
392     yu1 = 0.0
393     yv1 = 0.0
394     yrads = 0.0
395     ypaprs = 0.0
396     ypplay = 0.0
397     ydelp = 0.0
398     yu = 0.0
399     yv = 0.0
400     yt = 0.0
401     yq = 0.0
402     pctsrf_new = 0.0
403     y_flux_u = 0.0
404     y_flux_v = 0.0
405     C$$ PB
406     y_dflux_t = 0.0
407     y_dflux_q = 0.0
408     ytsoil = 999999.
409     yrugoro = 0.
410     c -- LOOP
411     yu10mx = 0.0
412     yu10my = 0.0
413     ywindsp = 0.0
414     c -- LOOP
415     DO nsrf = 1, nbsrf
416     DO i = 1, klon
417     d_ts(i,nsrf) = 0.0
418     ENDDO
419     END DO
420     C§§§ PB
421     yfluxlat=0.
422     flux_t = 0.
423     flux_q = 0.
424     flux_u = 0.
425     flux_v = 0.
426     DO k = 1, klev
427     DO i = 1, klon
428     d_t(i,k) = 0.0
429     d_q(i,k) = 0.0
430     c$$$ flux_t(i,k) = 0.0
431     c$$$ flux_q(i,k) = 0.0
432     d_u(i,k) = 0.0
433     d_v(i,k) = 0.0
434     c$$$ flux_u(i,k) = 0.0
435     c$$$ flux_v(i,k) = 0.0
436     zcoefh(i,k) = 0.0
437     ENDDO
438     ENDDO
439     cAA IF (itr.GE.1) THEN
440     cAA DO it = 1, itr
441     cAA DO k = 1, klev
442     cAA DO i = 1, klon
443     cAA d_tr(i,k,it) = 0.0
444     cAA ENDDO
445     cAA ENDDO
446     cAA ENDDO
447     cAA ENDIF
448    
449     c
450     c Boucler sur toutes les sous-fractions du sol:
451     c
452     C Initialisation des "pourcentages potentiels". On considere ici qu'on
453     C peut avoir potentiellementdela glace sur tout le domaine oceanique
454     C (a affiner)
455    
456     pctsrf_pot = pctsrf
457     pctsrf_pot(:,is_oce) = 1. - zmasq(:)
458     pctsrf_pot(:,is_sic) = 1. - zmasq(:)
459    
460     DO 99999 nsrf = 1, nbsrf
461    
462     c chercher les indices:
463     DO j = 1, klon
464     ni(j) = 0
465     ENDDO
466     knon = 0
467     DO i = 1, klon
468    
469     C pour determiner le domaine a traiter on utilise les surfaces "potentielles"
470     C
471     IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN
472     knon = knon + 1
473     ni(knon) = i
474     ENDIF
475     ENDDO
476     c
477     if (check) THEN
478     write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon
479     CC call flush(6)
480     endif
481     c
482     c variables pour avoir une sortie IOIPSL des INDEX
483     c
484     IF (debugindex) THEN
485     tabindx(:)=0.
486     c tabindx(1:knon)=(/FLOAT(i),i=1:knon/)
487     DO i=1,knon
488     tabindx(1:knon)=FLOAT(i)
489     END DO
490     debugtab(:,:)=0.
491     ndexbg(:)=0
492     CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)
493     CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)
494     $ ,ndexbg)
495     ENDIF
496     IF (knon.EQ.0) GOTO 99999
497     DO j = 1, knon
498     i = ni(j)
499     ypct(j) = pctsrf(i,nsrf)
500     yts(j) = ts(i,nsrf)
501     cIM "slab" ocean
502     c PRINT *, 'tslab = ', i, tslab(i)
503     ytslab(i) = tslab(i)
504     c
505     ysnow(j) = snow(i,nsrf)
506     yqsurf(j) = qsurf(i,nsrf)
507     yalb(j) = albe(i,nsrf)
508     yalblw(j) = alblw(i,nsrf)
509     yrain_f(j) = rain_f(i)
510     ysnow_f(j) = snow_f(i)
511     yagesno(j) = agesno(i,nsrf)
512     yfder(j) = fder(i)
513     ytaux(j) = flux_u(i,1,nsrf)
514     ytauy(j) = flux_v(i,1,nsrf)
515     ysolsw(j) = solsw(i,nsrf)
516     ysollw(j) = sollw(i,nsrf)
517     ysollwdown(j) = sollwdown(i)
518     yrugos(j) = rugos(i,nsrf)
519     yrugoro(j) = rugoro(i)
520     yu1(j) = u1lay(i)
521     yv1(j) = v1lay(i)
522     yrads(j) = ysolsw(j)+ ysollw(j)
523     ypaprs(j,klev+1) = paprs(i,klev+1)
524     y_run_off_lic_0(j) = run_off_lic_0(i)
525     c -- LOOP
526     yu10mx(j) = u10m(i,nsrf)
527     yu10my(j) = v10m(i,nsrf)
528     ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )
529     c -- LOOP
530     END DO
531     C
532     C IF bucket model for continent, copy soil water content
533     IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN
534     DO j = 1, knon
535     i = ni(j)
536     yqsol(j) = qsol(i)
537     END DO
538     ELSE
539     yqsol(:)=0.
540     ENDIF
541     c$$$ PB ajour pour soil
542     DO k = 1, nsoilmx
543     DO j = 1, knon
544     i = ni(j)
545     ytsoil(j,k) = ftsoil(i,k,nsrf)
546     END DO
547     END DO
548     DO k = 1, klev
549     DO j = 1, knon
550     i = ni(j)
551     ypaprs(j,k) = paprs(i,k)
552     ypplay(j,k) = pplay(i,k)
553     ydelp(j,k) = delp(i,k)
554     yu(j,k) = u(i,k)
555     yv(j,k) = v(i,k)
556     yt(j,k) = t(i,k)
557     yq(j,k) = q(i,k)
558     ENDDO
559     ENDDO
560     c
561     c
562     c calculer Cdrag et les coefficients d'echange
563     CALL coefkz(nsrf, knon, ypaprs, ypplay,
564     cIM 261103
565     . ksta, ksta_ter,
566     cIM 261103
567     . yts, yrugos, yu, yv, yt, yq,
568     . yqsurf,
569     . ycoefm, ycoefh)
570     cIM 081204 BEG
571     cCR test
572     if (iflag_pbl.eq.1) then
573     cIM 081204 END
574     CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,
575     . ycoefm0, ycoefh0)
576     DO k = 1, klev
577     DO i = 1, knon
578     ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
579     ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
580     ENDDO
581     ENDDO
582     endif
583     c
584     cIM cf JLD : on seuille ycoefm et ycoefh
585     if (nsrf.eq.is_oce) then
586     do j=1,knon
587     c ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)
588     ycoefm(j,1)=min(ycoefm(j,1),cdmmax)
589     c ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)
590     ycoefh(j,1)=min(ycoefh(j,1),cdhmax)
591     enddo
592     endif
593    
594     c
595     cIM: 261103
596     if (ok_kzmin) THEN
597     cIM cf FH: 201103 BEG
598     c Calcul d'une diffusion minimale pour les conditions tres stables.
599     call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm
600     . ,ycoefm0,ycoefh0)
601     c call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ ')
602     c call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN ')
603    
604     if ( 1.eq.1 ) then
605     DO k = 1, klev
606     DO i = 1, knon
607     ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
608     ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))
609     ENDDO
610     ENDDO
611     endif
612     cIM cf FH: 201103 END
613     endif !ok_kzmin
614     cIM: 261103
615    
616    
617     IF (iflag_pbl.ge.3) then
618    
619     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
620     c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin
621     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
622    
623     yzlay(1:knon,1)=
624     . RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1)))
625     . *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG
626     do k=2,klev
627     yzlay(1:knon,k)=
628     . yzlay(1:knon,k-1)+RD*0.5*(yt(1:knon,k-1)+yt(1:knon,k))
629     . /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG
630     enddo
631     do k=1,klev
632     yteta(1:knon,k)=
633     . yt(1:knon,k)*(ypaprs(1:knon,1)/ypplay(1:knon,k))**rkappa
634     . *(1.+0.61*yq(1:knon,k))
635     enddo
636     yzlev(1:knon,1)=0.
637     yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)
638     do k=2,klev
639     yzlev(1:knon,k)=0.5*(yzlay(1:knon,k)+yzlay(1:knon,k-1))
640     enddo
641     DO k = 1, klev+1
642     DO j = 1, knon
643     i = ni(j)
644     yq2(j,k)=q2(i,k,nsrf)
645     enddo
646     enddo
647    
648    
649     c Bug introduit volontairement pour converger avec les resultats
650     c du papier sur les thermiques.
651     if (1.eq.1) then
652     y_cd_m(1:knon) = ycoefm(1:knon,1)
653     y_cd_h(1:knon) = ycoefh(1:knon,1)
654     else
655     y_cd_h(1:knon) = ycoefm(1:knon,1)
656     y_cd_m(1:knon) = ycoefh(1:knon,1)
657     endif
658     call ustarhb(knon,yu,yv,y_cd_m, yustar)
659    
660     if (prt_level > 9) THEN
661     WRITE(lunout,*)'USTAR = ',yustar
662     ENDIF
663    
664     c iflag_pbl peut etre utilise comme longuer de melange
665    
666     if (iflag_pbl.ge.11) then
667     call vdif_kcay(knon,dtime,rg,rd,ypaprs,yt
668     s ,yzlev,yzlay,yu,yv,yteta
669     s ,y_cd_m,yq2,q2diag,ykmm,ykmn,yustar,
670     s iflag_pbl)
671     else
672     call yamada4(knon,dtime,rg,rd,ypaprs,yt
673     s ,yzlev,yzlay,yu,yv,yteta
674     s ,y_cd_m,yq2,ykmm,ykmn,ykmq,yustar,
675     s iflag_pbl)
676     endif
677    
678     ycoefm(1:knon,1)=y_cd_m(1:knon)
679     ycoefh(1:knon,1)=y_cd_h(1:knon)
680     ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
681     ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
682    
683    
684     ENDIF
685    
686     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
687     c calculer la diffusion des vitesses "u" et "v"
688     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
689    
690     CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
691     s y_d_u,y_flux_u)
692     CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
693     s y_d_v,y_flux_v)
694    
695     c pour le couplage
696     ytaux = y_flux_u(:,1)
697     ytauy = y_flux_v(:,1)
698    
699     c FH modif sur le cdrag temperature
700     c$$$PB : déplace dans clcdrag
701     c$$$ do i=1,knon
702     c$$$ ycoefh(i,1)=ycoefm(i,1)*0.8
703     c$$$ enddo
704    
705     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
706     c calculer la diffusion de "q" et de "h"
707     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
708     CALL clqh(dtime, itap, date0,jour, debut,lafin,
709     e rlon, rlat, cufi, cvfi,
710     e knon, nsrf, ni, pctsrf,
711     e soil_model, ytsoil,yqsol,
712     e ok_veget, ocean, npas, nexca,
713     e rmu0, co2_ppm, yrugos, yrugoro,
714     e yu1, yv1, ycoefh,
715     e yt,yq,yts,ypaprs,ypplay,
716     e ydelp,yrads,yalb, yalblw, ysnow, yqsurf,
717     e yrain_f, ysnow_f, yfder, ytaux, ytauy,
718     c -- LOOP
719     e ywindsp,
720     c -- LOOP
721     c$$$ e ysollw, ysolsw,
722     e ysollw, ysollwdown, ysolsw,yfluxlat,
723     s pctsrf_new, yagesno,
724     s y_d_t, y_d_q, y_d_ts, yz0_new,
725     s y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,
726     s y_fqcalving,y_ffonte,y_run_off_lic_0,
727     cIM "slab" ocean
728     s y_flux_o, y_flux_g, ytslab, y_seaice)
729     c
730     c calculer la longueur de rugosite sur ocean
731     yrugm=0.
732     IF (nsrf.EQ.is_oce) THEN
733     DO j = 1, knon
734     yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG
735     $ + 0.11*14e-6 / sqrt(ycoefm(j,1) * (yu1(j)**2+yv1(j)**2))
736     yrugm(j) = MAX(1.5e-05,yrugm(j))
737     ENDDO
738     ENDIF
739     DO j = 1, knon
740     y_dflux_t(j) = y_dflux_t(j) * ypct(j)
741     y_dflux_q(j) = y_dflux_q(j) * ypct(j)
742     yu1(j) = yu1(j) * ypct(j)
743     yv1(j) = yv1(j) * ypct(j)
744     ENDDO
745     c
746     DO k = 1, klev
747     DO j = 1, knon
748     i = ni(j)
749     ycoefh(j,k) = ycoefh(j,k) * ypct(j)
750     ycoefm(j,k) = ycoefm(j,k) * ypct(j)
751     y_d_t(j,k) = y_d_t(j,k) * ypct(j)
752     y_d_q(j,k) = y_d_q(j,k) * ypct(j)
753     C§§§ PB
754     flux_t(i,k,nsrf) = y_flux_t(j,k)
755     flux_q(i,k,nsrf) = y_flux_q(j,k)
756     flux_u(i,k,nsrf) = y_flux_u(j,k)
757     flux_v(i,k,nsrf) = y_flux_v(j,k)
758     c$$$ PB y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)
759     c$$$ PB y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)
760     y_d_u(j,k) = y_d_u(j,k) * ypct(j)
761     y_d_v(j,k) = y_d_v(j,k) * ypct(j)
762     c$$$ PB y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)
763     c$$$ PB y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)
764     ENDDO
765     ENDDO
766    
767    
768     evap(:,nsrf) = - flux_q(:,1,nsrf)
769     c
770     albe(:, nsrf) = 0.
771     alblw(:, nsrf) = 0.
772     snow(:, nsrf) = 0.
773     qsurf(:, nsrf) = 0.
774     rugos(:, nsrf) = 0.
775     fluxlat(:,nsrf) = 0.
776     DO j = 1, knon
777     i = ni(j)
778     d_ts(i,nsrf) = y_d_ts(j)
779     albe(i,nsrf) = yalb(j)
780     alblw(i,nsrf) = yalblw(j)
781     snow(i,nsrf) = ysnow(j)
782     qsurf(i,nsrf) = yqsurf(j)
783     rugos(i,nsrf) = yz0_new(j)
784     fluxlat(i,nsrf) = yfluxlat(j)
785     c$$$ pb rugmer(i) = yrugm(j)
786     IF (nsrf .EQ. is_oce) then
787     rugmer(i) = yrugm(j)
788     rugos(i,nsrf) = yrugm(j)
789     endif
790     cIM cf JLD ??
791     agesno(i,nsrf) = yagesno(j)
792     fqcalving(i,nsrf) = y_fqcalving(j)
793     ffonte(i,nsrf) = y_ffonte(j)
794     cdragh(i) = cdragh(i) + ycoefh(j,1)
795     cdragm(i) = cdragm(i) + ycoefm(j,1)
796     dflux_t(i) = dflux_t(i) + y_dflux_t(j)
797     dflux_q(i) = dflux_q(i) + y_dflux_q(j)
798     zu1(i) = zu1(i) + yu1(j)
799     zv1(i) = zv1(i) + yv1(j)
800     END DO
801     IF ( nsrf .eq. is_ter ) THEN
802     DO j = 1, knon
803     i = ni(j)
804     qsol(i) = yqsol(j)
805     END DO
806     END IF
807     IF ( nsrf .eq. is_lic ) THEN
808     DO j = 1, knon
809     i = ni(j)
810     run_off_lic_0(i) = y_run_off_lic_0(j)
811     END DO
812     END IF
813     c$$$ PB ajout pour soil
814     ftsoil(:,:,nsrf) = 0.
815     DO k = 1, nsoilmx
816     DO j = 1, knon
817     i = ni(j)
818     ftsoil(i, k, nsrf) = ytsoil(j,k)
819     END DO
820     END DO
821     c
822     DO j = 1, knon
823     i = ni(j)
824     DO k = 1, klev
825     d_t(i,k) = d_t(i,k) + y_d_t(j,k)
826     d_q(i,k) = d_q(i,k) + y_d_q(j,k)
827     c$$$ PB flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)
828     c$$$ flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)
829     d_u(i,k) = d_u(i,k) + y_d_u(j,k)
830     d_v(i,k) = d_v(i,k) + y_d_v(j,k)
831     c$$$ PB flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)
832     c$$$ flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)
833     zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
834     ENDDO
835     ENDDO
836     c
837     c
838     ccc diagnostic t,q a 2m et u, v a 10m
839     c
840     DO j=1, knon
841     i = ni(j)
842     uzon(j) = yu(j,1) + y_d_u(j,1)
843     vmer(j) = yv(j,1) + y_d_v(j,1)
844     tair1(j) = yt(j,1) + y_d_t(j,1)
845     qair1(j) = yq(j,1) + y_d_q(j,1)
846     zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1)))
847     & * (ypaprs(j,1)-ypplay(j,1))
848     tairsol(j) = yts(j) + y_d_ts(j)
849     rugo1(j) = yrugos(j)
850     IF(nsrf.EQ.is_oce) THEN
851     rugo1(j) = rugos(i,nsrf)
852     ENDIF
853     psfce(j)=ypaprs(j,1)
854     patm(j)=ypplay(j,1)
855     c
856     qairsol(j) = yqsurf(j)
857     ENDDO
858     c
859     CALL stdlevvar(klon, knon, nsrf, zxli,
860     & uzon, vmer, tair1, qair1, zgeo1,
861     & tairsol, qairsol, rugo1, psfce, patm,
862     cIM & yt2m, yq2m, yu10m)
863     & yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
864     cIM 081204 END
865     c
866     c
867     DO j=1, knon
868     i = ni(j)
869     t2m(i,nsrf)=yt2m(j)
870    
871     c
872     q2m(i,nsrf)=yq2m(j)
873     c
874     c u10m, v10m : composantes du vent a 10m sans spirale de Ekman
875     u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
876     v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
877     c
878     ENDDO
879     c
880     cIM cf AM : pbl, HBTM
881     DO i = 1, knon
882     y_cd_h(i) = ycoefh(i,1)
883     y_cd_m(i) = ycoefm(i,1)
884     ENDDO
885     c print*,'appel hbtm2'
886     CALL HBTM(knon, ypaprs, ypplay,
887     . yt2m,yt10m,yq2m,yq10m,yustar,
888     . y_flux_t,y_flux_q,yu,yv,yt,yq,
889     . ypblh,ycapCL,yoliqCL,ycteiCL,ypblT,
890     . ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
891     c print*,'fin hbtm2'
892     c
893     DO j=1, knon
894     i = ni(j)
895     pblh(i,nsrf) = ypblh(j)
896     plcl(i,nsrf) = ylcl(j)
897     capCL(i,nsrf) = ycapCL(j)
898     oliqCL(i,nsrf) = yoliqCL(j)
899     cteiCL(i,nsrf) = ycteiCL(j)
900     pblT(i,nsrf) = ypblT(j)
901     therm(i,nsrf) = ytherm(j)
902     trmb1(i,nsrf) = ytrmb1(j)
903     trmb2(i,nsrf) = ytrmb2(j)
904     trmb3(i,nsrf) = ytrmb3(j)
905     ENDDO
906     c
907    
908     do j=1,knon
909     do k=1,klev+1
910     i=ni(j)
911     q2(i,k,nsrf)=yq2(j,k)
912     enddo
913     enddo
914     cIM "slab" ocean
915     IF(OCEAN.EQ.'slab '.OR.OCEAN.EQ.'force ') THEN
916     IF (nsrf.EQ.is_oce) THEN
917     DO j = 1, knon
918     c on projette sur la grille globale
919     i = ni(j)
920     IF(pctsrf_new(i,is_oce).GT.epsfra) THEN
921     flux_o(i) = y_flux_o(j)
922     ELSE
923     flux_o(i) = 0.
924     ENDIF
925     ENDDO
926     ENDIF
927     c
928     IF (nsrf.EQ.is_sic) THEN
929     DO j = 1, knon
930     i = ni(j)
931     cIM 230604 on pondere lorsque l'on fait le bilan au sol : flux_g(i) = y_flux_g(j)*ypct(j)
932     IF(pctsrf_new(i,is_sic).GT.epsfra) THEN
933     flux_g(i) = y_flux_g(j)
934     ELSE
935     flux_g(i) = 0.
936     ENDIF
937     ENDDO
938     ENDIF !nsrf.EQ.is_sic
939     ENDIF !OCEAN
940     c
941     IF(OCEAN.EQ.'slab ') THEN
942     IF(nsrf.EQ.is_oce) then
943     tslab(1:klon) = ytslab(1:klon)
944     seaice(1:klon) = y_seaice(1:klon)
945     ENDIF !nsrf
946     ENDIF !OCEAN
947     99999 CONTINUE
948     C
949     C On utilise les nouvelles surfaces
950     C A rajouter: conservation de l'albedo
951     C
952     rugos(:,is_oce) = rugmer
953     pctsrf = pctsrf_new
954    
955     RETURN
956     END

  ViewVC Help
Powered by ViewVC 1.1.21