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

Diff of /trunk/phylmd/clmain.f

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

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

Legend:
Removed from v.3  
changed lines
  Added in v.30

  ViewVC Help
Powered by ViewVC 1.1.21