/[lmdze]/trunk/Sources/phylmd/clmain.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/clmain.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

  ViewVC Help
Powered by ViewVC 1.1.21