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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21