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

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

  ViewVC Help
Powered by ViewVC 1.1.21