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

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

  ViewVC Help
Powered by ViewVC 1.1.21