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

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

  ViewVC Help
Powered by ViewVC 1.1.21