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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 10 - (hide 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 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 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     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