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

Legend:
Removed from v.10  
changed lines
  Added in v.37

  ViewVC Help
Powered by ViewVC 1.1.21