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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21