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

Legend:
Removed from v.13  
changed lines
  Added in v.61

  ViewVC Help
Powered by ViewVC 1.1.21