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

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

  ViewVC Help
Powered by ViewVC 1.1.21