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

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

  ViewVC Help
Powered by ViewVC 1.1.21