/[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 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, intent(in):: 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)    REAL d_ts(klon, nbsrf)
145        real, intent(in):: rugoro(klon)    REAL snow(klon, nbsrf)
146        REAL cdragh(klon), cdragm(klon)    REAL qsurf(klon, nbsrf)
147        integer jour            ! jour de l'annee en cours    REAL evap(klon, nbsrf)
148        real rmu0(klon)         ! cosinus de l'angle solaire zenithal    REAL albe(klon, nbsrf)
149        REAL co2_ppm            ! taux CO2 atmosphere    REAL alblw(klon, nbsrf)
150        LOGICAL, intent(in):: debut  
151        logical, intent(in):: lafin    REAL fluxlat(klon, nbsrf)
152        logical ok_veget  
153        character(len=*), intent(IN):: ocean    REAL rain_f(klon), snow_f(klon)
154        integer npas, nexca    REAL fder(klon)
155  c  
156        REAL pctsrf(klon,nbsrf)    REAL sollw(klon, nbsrf), solsw(klon, nbsrf), sollwdown(klon)
157        REAL ts(klon,nbsrf)    REAL rugos(klon, nbsrf)
158        REAL d_ts(klon,nbsrf)    ! la nouvelle repartition des surfaces sortie de l'interface
159        REAL snow(klon,nbsrf)    REAL pctsrf_new(klon, nbsrf)
160        REAL qsurf(klon,nbsrf)  
161        REAL evap(klon,nbsrf)    REAL zcoefh(klon, klev)
162        REAL albe(klon,nbsrf)    REAL zu1(klon)
163        REAL alblw(klon,nbsrf)    REAL zv1(klon)
164  c$$$ PB  
165        REAL fluxlat(klon,nbsrf)    !$$$ PB ajout pour soil
166  C    LOGICAL, INTENT (IN) :: soil_model
167        real rain_f(klon), snow_f(klon)    !IM ajout seuils cdrm, cdrh
168        REAL fder(klon)    REAL cdmmax, cdhmax
169  cIM cf. JLD   REAL sollw(klon), solsw(klon), sollwdown(klon)  
170        REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)    REAL ksta, ksta_ter
171        REAL rugos(klon,nbsrf)    LOGICAL ok_kzmin
172  C la nouvelle repartition des surfaces sortie de l'interface  
173        REAL pctsrf_new(klon,nbsrf)    REAL ftsoil(klon, nsoilmx, nbsrf)
174  cAA    REAL ytsoil(klon, nsoilmx)
175        REAL zcoefh(klon,klev)    REAL qsol(klon)
176        REAL zu1(klon)  
177        REAL zv1(klon)    EXTERNAL clqh, clvent, coefkz, calbeta, cltrac
178  cAA  
179  c$$$ PB ajout pour soil    REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)
180        LOGICAL, intent(in):: soil_model    REAL yalb(klon)
181  cIM ajout seuils cdrm, cdrh    REAL yalblw(klon)
182        REAL cdmmax, cdhmax    REAL yu1(klon), yv1(klon)
183  cIM: 261103    REAL ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)
184        REAL ksta, ksta_ter    REAL yrain_f(klon), ysnow_f(klon)
185        LOGICAL ok_kzmin    REAL ysollw(klon), ysolsw(klon), ysollwdown(klon)
186  cIM: 261103    REAL yfder(klon), ytaux(klon), ytauy(klon)
187        REAL ftsoil(klon,nsoilmx,nbsrf)    REAL yrugm(klon), yrads(klon), yrugoro(klon)
188        REAL ytsoil(klon,nsoilmx)  
189        REAL qsol(klon)    REAL yfluxlat(klon)
190  c======================================================================  
191        EXTERNAL clqh, clvent, coefkz, calbeta, cltrac    REAL y_d_ts(klon)
192  c======================================================================    REAL y_d_t(klon, klev), y_d_q(klon, klev)
193        REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)    REAL y_d_u(klon, klev), y_d_v(klon, klev)
194        REAL yalb(klon)    REAL y_flux_t(klon, klev), y_flux_q(klon, klev)
195        REAL yalblw(klon)    REAL y_flux_u(klon, klev), y_flux_v(klon, klev)
196        REAL yu1(klon), yv1(klon)    REAL y_dflux_t(klon), y_dflux_q(klon)
197        real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)    REAL ycoefh(klon, klev), ycoefm(klon, klev)
198        real yrain_f(klon), ysnow_f(klon)    REAL yu(klon, klev), yv(klon, klev)
199        real ysollw(klon), ysolsw(klon), ysollwdown(klon)    REAL yt(klon, klev), yq(klon, klev)
200        real yfder(klon), ytaux(klon), ytauy(klon)    REAL ypaprs(klon, klev+1), ypplay(klon, klev), ydelp(klon, klev)
201        REAL yrugm(klon), yrads(klon),yrugoro(klon)  
202  c$$$ PB    LOGICAL ok_nonloc
203        REAL yfluxlat(klon)    PARAMETER (ok_nonloc=.FALSE.)
204  C    REAL ycoefm0(klon, klev), ycoefh0(klon, klev)
205        REAL y_d_ts(klon)  
206        REAL y_d_t(klon, klev), y_d_q(klon, klev)    !IM 081204 hcl_Anne ? BEG
207        REAL y_d_u(klon, klev), y_d_v(klon, klev)    REAL yzlay(klon, klev), yzlev(klon, klev+1), yteta(klon, klev)
208        REAL y_flux_t(klon,klev), y_flux_q(klon,klev)    REAL ykmm(klon, klev+1), ykmn(klon, klev+1)
209        REAL y_flux_u(klon,klev), y_flux_v(klon,klev)    REAL ykmq(klon, klev+1)
210        REAL y_dflux_t(klon), y_dflux_q(klon)    REAL yq2(klon, klev+1), q2(klon, klev+1, nbsrf)
211        REAL ycoefh(klon,klev), ycoefm(klon,klev)    REAL q2diag(klon, klev+1)
212        REAL yu(klon,klev), yv(klon,klev)    !IM 081204 hcl_Anne ? END
213        REAL yt(klon,klev), yq(klon,klev)  
214        REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)    REAL u1lay(klon), v1lay(klon)
215  c    REAL delp(klon, klev)
216        LOGICAL ok_nonloc    INTEGER i, k, nsrf
217        PARAMETER (ok_nonloc=.FALSE.)  
218        REAL ycoefm0(klon,klev), ycoefh0(klon,klev)    INTEGER ni(klon), knon, j
219      ! Introduction d'une variable "pourcentage potentiel" pour tenir compte
220  cIM 081204 hcl_Anne ? BEG    ! des eventuelles apparitions et/ou disparitions de la glace de mer
221        real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)    REAL pctsrf_pot(klon, nbsrf)
222        real ykmm(klon,klev+1),ykmn(klon,klev+1)  
223        real ykmq(klon,klev+1)    REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.
224        real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)  
225        real q2diag(klon,klev+1)    ! maf pour sorties IOISPL en cas de debugagage
226  cIM 081204   real yustar(klon),y_cd_m(klon),y_cd_h(klon)  
227  cIM 081204 hcl_Anne ? END    CHARACTER (80) cldebug
228  c    SAVE cldebug
229        REAL u1lay(klon), v1lay(klon)    CHARACTER (8) cl_surf(nbsrf)
230        REAL delp(klon,klev)    SAVE cl_surf
231        INTEGER i, k, nsrf    INTEGER nhoridbg, nidbg
232  cAA   INTEGER it    SAVE nhoridbg, nidbg
233        INTEGER ni(klon), knon, j    INTEGER ndexbg(iim*(jjm+1))
234  c Introduction d'une variable "pourcentage potentiel" pour tenir compte    REAL zx_lon(iim, jjm+1), zx_lat(iim, jjm+1), zjulian
235  c des eventuelles apparitions et/ou disparitions de la glace de mer    REAL tabindx(klon)
236        REAL pctsrf_pot(klon,nbsrf)    REAL debugtab(iim, jjm+1)
237            LOGICAL first_appel
238  c======================================================================    SAVE first_appel
239        REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.    DATA first_appel/ .TRUE./
240  c======================================================================    LOGICAL :: debugindex = .FALSE.
241  c    INTEGER idayref
242  c maf pour sorties IOISPL en cas de debugagage    REAL t2m(klon, nbsrf), q2m(klon, nbsrf)
243  c    REAL u10m(klon, nbsrf), v10m(klon, nbsrf)
244        CHARACTER*80 cldebug  
245        SAVE cldebug    REAL yt2m(klon), yq2m(klon), yu10m(klon)
246        CHARACTER*8 cl_surf(nbsrf)    REAL yustar(klon)
247        SAVE cl_surf    ! -- LOOP
248        INTEGER nhoridbg, nidbg    REAL yu10mx(klon)
249        SAVE nhoridbg, nidbg    REAL yu10my(klon)
250        INTEGER ndexbg(iim*(jjm+1))    REAL ywindsp(klon)
251        REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian    ! -- LOOP
252        REAL tabindx(klon)  
253        REAL debugtab(iim,jjm+1)    REAL yt10m(klon), yq10m(klon)
254        LOGICAL first_appel    !IM cf. AM : pbl, hbtm (Comme les autres diagnostics on cumule ds
255        SAVE first_appel    ! physiq ce qui permet de sortir les grdeurs par sous surface)
256        DATA first_appel/.true./    REAL pblh(klon, nbsrf)
257        LOGICAL:: debugindex = .false.    REAL plcl(klon, nbsrf)
258        integer idayref    REAL capcl(klon, nbsrf)
259        REAL t2m(klon,nbsrf), q2m(klon,nbsrf)    REAL oliqcl(klon, nbsrf)
260        REAL u10m(klon,nbsrf), v10m(klon,nbsrf)    REAL cteicl(klon, nbsrf)
261  c    REAL pblt(klon, nbsrf)
262        REAL yt2m(klon), yq2m(klon), yu10m(klon)    REAL therm(klon, nbsrf)
263        REAL yustar(klon)    REAL trmb1(klon, nbsrf)
264  c -- LOOP    REAL trmb2(klon, nbsrf)
265         REAL yu10mx(klon)    REAL trmb3(klon, nbsrf)
266         REAL yu10my(klon)    REAL ypblh(klon)
267         REAL ywindsp(klon)    REAL ylcl(klon)
268  c -- LOOP    REAL ycapcl(klon)
269  c    REAL yoliqcl(klon)
270        REAL yt10m(klon), yq10m(klon)    REAL ycteicl(klon)
271  cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui    REAL ypblt(klon)
272  c   permet de sortir les grdeurs par sous surface)    REAL ytherm(klon)
273        REAL pblh(klon,nbsrf)    REAL ytrmb1(klon)
274        REAL plcl(klon,nbsrf)    REAL ytrmb2(klon)
275        REAL capCL(klon,nbsrf)    REAL ytrmb3(klon)
276        REAL oliqCL(klon,nbsrf)    REAL y_cd_h(klon), y_cd_m(klon)
277        REAL cteiCL(klon,nbsrf)    REAL uzon(klon), vmer(klon)
278        REAL pblT(klon,nbsrf)    REAL tair1(klon), qair1(klon), tairsol(klon)
279        REAL therm(klon,nbsrf)    REAL psfce(klon), patm(klon)
280        REAL trmb1(klon,nbsrf)  
281        REAL trmb2(klon,nbsrf)    REAL qairsol(klon), zgeo1(klon)
282        REAL trmb3(klon,nbsrf)    REAL rugo1(klon)
283        REAL ypblh(klon)  
284        REAL ylcl(klon)    ! utiliser un jeu de fonctions simples              
285        REAL ycapCL(klon)    LOGICAL zxli
286        REAL yoliqCL(klon)    PARAMETER (zxli=.FALSE.)
287        REAL ycteiCL(klon)  
288        REAL ypblT(klon)    REAL zt, zqs, zdelta, zcor
289        REAL ytherm(klon)    REAL t_coup
290        REAL ytrmb1(klon)    PARAMETER (t_coup=273.15)
291        REAL ytrmb2(klon)  
292        REAL ytrmb3(klon)    CHARACTER (len=20) :: modname = 'clmain'
293        REAL y_cd_h(klon), y_cd_m(klon)  
294  c     REAL ygamt(klon,2:klev) ! contre-gradient pour temperature    !------------------------------------------------------------
295  c     REAL ygamq(klon,2:klev) ! contre-gradient pour humidite  
296        REAL uzon(klon), vmer(klon)    ! initialisation Anne
297        REAL tair1(klon), qair1(klon), tairsol(klon)    ytherm = 0.
298        REAL psfce(klon), patm(klon)  
299  c    IF (debugindex .AND. first_appel) THEN
300        REAL qairsol(klon), zgeo1(klon)       first_appel = .FALSE.
301        REAL rugo1(klon)  
302  c       ! initialisation sorties netcdf
303        LOGICAL zxli ! utiliser un jeu de fonctions simples  
304        PARAMETER (zxli=.FALSE.)       idayref = day_ini
305  c       CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
306        REAL zt, zqs, zdelta, zcor       CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlon, zx_lon)
307        REAL t_coup       DO i = 1, iim
308        PARAMETER(t_coup=273.15)          zx_lon(i, 1) = rlon(i+1)
309  C          zx_lon(i, jjm+1) = rlon(i+1)
310        character (len = 20) :: modname = 'clmain'       END DO
311        LOGICAL check       CALL gr_fi_ecrit(1, klon, iim, jjm+1, rlat, zx_lat)
312        PARAMETER (check=.false.)       cldebug = 'sous_index'
313         CALL histbeg_totreg(cldebug, zx_lon(:, 1), zx_lat(1, :), 1, &
314              iim, 1, jjm+1, itau_phy, zjulian, dtime, nhoridbg, nidbg)
315  c initialisation Anne       ! no vertical axis
316        ytherm(:) = 0.       cl_surf(1) = 'ter'
317  C       cl_surf(2) = 'lic'
318        if (check) THEN       cl_surf(3) = 'oce'
319            write(*,*) modname,'  klon=',klon       cl_surf(4) = 'sic'
320  CC        call flush(6)       DO nsrf = 1, nbsrf
321        endif          CALL histdef(nidbg, cl_surf(nsrf), cl_surf(nsrf), '-', iim, jjm+1, &
322        IF (debugindex .and. first_appel) THEN               nhoridbg, 1, 1, 1, -99, 'inst', dtime, dtime)
323            first_appel=.false.       END DO
324  !       CALL histend(nidbg)
325  ! initialisation sorties netcdf       CALL histsync(nidbg)
326  !    END IF
327            idayref = day_ini  
328            CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)    DO k = 1, klev ! epaisseur de couche
329            CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)       DO i = 1, klon
330            DO i = 1, iim          delp(i, k) = paprs(i, k) - paprs(i, k+1)
331              zx_lon(i,1) = rlon(i+1)       END DO
332              zx_lon(i,jjm+1) = rlon(i+1)    END DO
333            ENDDO    DO i = 1, klon ! vent de la premiere couche
334            CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)       zx_alf1 = 1.0
335            cldebug='sous_index'       zx_alf2 = 1.0 - zx_alf1
336            CALL histbeg_totreg(cldebug, iim,zx_lon(:,1),jjm+1,       u1lay(i) = u(i, 1)*zx_alf1 + u(i, 2)*zx_alf2
337       $        zx_lat(1,:),1,iim,1,jjm       v1lay(i) = v(i, 1)*zx_alf1 + v(i, 2)*zx_alf2
338       $        +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)    END DO
339  ! no vertical axis  
340            cl_surf(1)='ter'    ! initialisation:
341            cl_surf(2)='lic'  
342            cl_surf(3)='oce'    DO i = 1, klon
343            cl_surf(4)='sic'       rugmer(i) = 0.0
344            DO nsrf=1,nbsrf       cdragh(i) = 0.0
345              CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,       cdragm(i) = 0.0
346       $          jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)       dflux_t(i) = 0.0
347            END DO       dflux_q(i) = 0.0
348            CALL histend(nidbg)       zu1(i) = 0.0
349            CALL histsync(nidbg)       zv1(i) = 0.0
350        ENDIF    END DO
351                ypct = 0.0
352        DO k = 1, klev   ! epaisseur de couche    yts = 0.0
353        DO i = 1, klon    ysnow = 0.0
354           delp(i,k) = paprs(i,k)-paprs(i,k+1)    yqsurf = 0.0
355        ENDDO    yalb = 0.0
356        ENDDO    yalblw = 0.0
357        DO i = 1, klon  ! vent de la premiere couche    yrain_f = 0.0
358           zx_alf1 = 1.0    ysnow_f = 0.0
359           zx_alf2 = 1.0 - zx_alf1    yfder = 0.0
360           u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2    ytaux = 0.0
361           v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2    ytauy = 0.0
362        ENDDO    ysolsw = 0.0
363  c    ysollw = 0.0
364  c initialisation:    ysollwdown = 0.0
365  c    yrugos = 0.0
366        DO i = 1, klon    yu1 = 0.0
367           rugmer(i) = 0.0    yv1 = 0.0
368           cdragh(i) = 0.0    yrads = 0.0
369           cdragm(i) = 0.0    ypaprs = 0.0
370           dflux_t(i) = 0.0    ypplay = 0.0
371           dflux_q(i) = 0.0    ydelp = 0.0
372           zu1(i) = 0.0    yu = 0.0
373           zv1(i) = 0.0    yv = 0.0
374        ENDDO    yt = 0.0
375        ypct = 0.0    yq = 0.0
376        yts = 0.0    pctsrf_new = 0.0
377        ysnow = 0.0    y_flux_u = 0.0
378        yqsurf = 0.0    y_flux_v = 0.0
379        yalb = 0.0    !$$ PB
380        yalblw = 0.0    y_dflux_t = 0.0
381        yrain_f = 0.0    y_dflux_q = 0.0
382        ysnow_f = 0.0    ytsoil = 999999.
383        yfder = 0.0    yrugoro = 0.
384        ytaux = 0.0    ! -- LOOP
385        ytauy = 0.0    yu10mx = 0.0
386        ysolsw = 0.0    yu10my = 0.0
387        ysollw = 0.0    ywindsp = 0.0
388        ysollwdown = 0.0    ! -- LOOP
389        yrugos = 0.0    DO nsrf = 1, nbsrf
390        yu1 = 0.0       DO i = 1, klon
391        yv1 = 0.0          d_ts(i, nsrf) = 0.0
392        yrads = 0.0       END DO
393        ypaprs = 0.0    END DO
394        ypplay = 0.0    !§§§ PB
395        ydelp = 0.0    yfluxlat = 0.
396        yu = 0.0    flux_t = 0.
397        yv = 0.0    flux_q = 0.
398        yt = 0.0    flux_u = 0.
399        yq = 0.0    flux_v = 0.
400        pctsrf_new = 0.0    DO k = 1, klev
401        y_flux_u = 0.0       DO i = 1, klon
402        y_flux_v = 0.0          d_t(i, k) = 0.0
403  C$$ PB          d_q(i, k) = 0.0
404        y_dflux_t = 0.0          !$$$         flux_t(i, k) = 0.0
405        y_dflux_q = 0.0          !$$$         flux_q(i, k) = 0.0
406        ytsoil = 999999.          d_u(i, k) = 0.0
407        yrugoro = 0.          d_v(i, k) = 0.0
408  c -- LOOP          !$$$         flux_u(i, k) = 0.0
409        yu10mx = 0.0          !$$$         flux_v(i, k) = 0.0
410        yu10my = 0.0          zcoefh(i, k) = 0.0
411        ywindsp = 0.0       END DO
412  c -- LOOP    END DO
413        DO nsrf = 1, nbsrf    !AA      IF (itr.GE.1) THEN
414        DO i = 1, klon    !AA      DO it = 1, itr
415           d_ts(i,nsrf) = 0.0    !AA      DO k = 1, klev
416        ENDDO    !AA      DO i = 1, klon
417        END DO    !AA         d_tr(i, k, it) = 0.0
418  C§§§ PB    !AA      ENDDO
419        yfluxlat=0.    !AA      ENDDO
420        flux_t = 0.    !AA      ENDDO
421        flux_q = 0.    !AA      ENDIF
422        flux_u = 0.  
423        flux_v = 0.  
424        DO k = 1, klev    ! Boucler sur toutes les sous-fractions du sol:
425        DO i = 1, klon  
426           d_t(i,k) = 0.0    ! Initialisation des "pourcentages potentiels". On considere ici qu'on
427           d_q(i,k) = 0.0    ! peut avoir potentiellementdela glace sur tout le domaine oceanique
428  c$$$         flux_t(i,k) = 0.0    ! (a affiner)
429  c$$$         flux_q(i,k) = 0.0  
430           d_u(i,k) = 0.0    pctsrf_pot = pctsrf
431           d_v(i,k) = 0.0    pctsrf_pot(:, is_oce) = 1. - zmasq
432  c$$$         flux_u(i,k) = 0.0    pctsrf_pot(:, is_sic) = 1. - zmasq
433  c$$$         flux_v(i,k) = 0.0  
434           zcoefh(i,k) = 0.0    DO nsrf = 1, nbsrf
435        ENDDO       ! chercher les indices:
436        ENDDO       ni = 0
437  cAA      IF (itr.GE.1) THEN       knon = 0
438  cAA      DO it = 1, itr       DO i = 1, klon
439  cAA      DO k = 1, klev          ! pour determiner le domaine a traiter on utilise les surfaces
440  cAA      DO i = 1, klon          ! "potentielles"
441  cAA         d_tr(i,k,it) = 0.0          IF (pctsrf_pot(i, nsrf) > epsfra) THEN
442  cAA      ENDDO             knon = knon + 1
443  cAA      ENDDO             ni(knon) = i
444  cAA      ENDDO          END IF
445  cAA      ENDIF       END DO
446    
447  c       ! variables pour avoir une sortie IOIPSL des INDEX
448  c Boucler sur toutes les sous-fractions du sol:       IF (debugindex) THEN
449  c          tabindx = 0.
450  C Initialisation des "pourcentages potentiels". On considere ici qu'on          DO i = 1, knon
451  C peut avoir potentiellementdela glace sur tout le domaine oceanique             tabindx(i) = real(i)
452  C (a affiner)          END DO
453            debugtab = 0.
454        pctsrf_pot = pctsrf          ndexbg = 0
455        pctsrf_pot(:,is_oce) = 1. - zmasq(:)          CALL gath2cpl(tabindx, debugtab, klon, knon, iim, jjm, ni)
456        pctsrf_pot(:,is_sic) = 1. - zmasq(:)          CALL histwrite(nidbg, cl_surf(nsrf), itap, debugtab)
457         END IF
458        DO 99999 nsrf = 1, nbsrf  
459         IF (knon==0) CYCLE
460  c chercher les indices:  
461        DO j = 1, klon       DO j = 1, knon
462           ni(j) = 0          i = ni(j)
463        ENDDO          ypct(j) = pctsrf(i, nsrf)
464        knon = 0          yts(j) = ts(i, nsrf)
       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            print *,'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 (nsrf.EQ.is_oce) THEN       END IF
849          DO j = 1, knon  
850  c on projette sur la grille globale       IF (nsrf==is_sic) THEN
851           i = ni(j)          DO j = 1, knon
852           IF(pctsrf_new(i,is_oce).GT.epsfra) THEN             i = ni(j)
853            flux_o(i) = y_flux_o(j)             !IM 230604 on pondere lorsque l'on fait le bilan au sol :  flux_g(i) = y_flux_g(j)*ypct(j)
854           ELSE             IF (pctsrf_new(i, is_sic)>epsfra) THEN
855            flux_o(i) = 0.                flux_g(i) = y_flux_g(j)
856           ENDIF             ELSE
857          ENDDO                flux_g(i) = 0.
858         ENDIF             END IF
859  c          END DO
860         IF (nsrf.EQ.is_sic) THEN  
861          DO j = 1, knon       END IF
862           i = ni(j)       !nsrf.EQ.is_sic                                            
863  cIM 230604 on pondere lorsque l'on fait le bilan au sol :  flux_g(i) = y_flux_g(j)*ypct(j)       IF (ocean=='slab  ') THEN
864           IF(pctsrf_new(i,is_sic).GT.epsfra) THEN          IF (nsrf==is_oce) THEN
865            flux_g(i) = y_flux_g(j)             tslab(1:klon) = ytslab(1:klon)
866           ELSE             seaice(1:klon) = y_seaice(1:klon)
867            flux_g(i) = 0.             !nsrf                                                      
868           ENDIF          END IF
869          ENDDO          !OCEAN                                                      
870         ENDIF !nsrf.EQ.is_sic       END IF
871  c    END DO
872        IF(OCEAN.EQ.'slab  ') THEN  
873         IF(nsrf.EQ.is_oce) then    ! On utilise les nouvelles surfaces
874          tslab(1:klon) = ytslab(1:klon)    ! A rajouter: conservation de l'albedo
875          seaice(1:klon) = y_seaice(1:klon)  
876         ENDIF !nsrf    rugos(:, is_oce) = rugmer
877        ENDIF !OCEAN    pctsrf = pctsrf_new
 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.13  
changed lines
  Added in v.37

  ViewVC Help
Powered by ViewVC 1.1.21