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

Legend:
Removed from v.14  
changed lines
  Added in v.15

  ViewVC Help
Powered by ViewVC 1.1.21