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

Legend:
Removed from v.12  
changed lines
  Added in v.49

  ViewVC Help
Powered by ViewVC 1.1.21