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

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

  ViewVC Help
Powered by ViewVC 1.1.21