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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.10  
changed lines
  Added in v.71

  ViewVC Help
Powered by ViewVC 1.1.21