/[lmdze]/trunk/libf/phylmd/clmain.f90
ViewVC logotype

Diff of /trunk/libf/phylmd/clmain.f90

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

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

Legend:
Removed from v.7  
changed lines
  Added in v.47

  ViewVC Help
Powered by ViewVC 1.1.21