/[lmdze]/trunk/phylmd/Interface_surf/pbl_surface.f90
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/pbl_surface.f90

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

trunk/libf/phylmd/clmain.f revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/phylmd/Interface_surf/pbl_surface.f90 revision 332 by guez, Tue Aug 13 09:19:22 2019 UTC
# Line 1  Line 1 
1        SUBROUTINE clmain(dtime,itap,date0,pctsrf,pctsrf_new,  module pbl_surface_m
2       .                  t,q,u,v,  
3       .                  jour, rmu0, co2_ppm,    IMPLICIT NONE
4       .                  ok_veget, ocean, npas, nexca, ts,  
5       .                  soil_model,cdmmax, cdhmax,  contains
6       .                  ksta, ksta_ter, ok_kzmin, ftsoil,qsol,  
7       .                  paprs,pplay,snow,qsurf,evap,albe,alblw,    SUBROUTINE pbl_surface(pctsrf, t, q, u, v, julien, mu0, ftsol, cdmmax, &
8       .                  fluxlat,         cdhmax, ftsoil, qsol, paprs, play, fsnow, fqsurf, falbe, fluxlat, &
9       .                  rain_f, snow_f, solsw, sollw, sollwdown, fder,         rain_fall, snow_fall, frugs, agesno, rugoro, d_t, d_q, d_u, d_v, &
10       .                  rlon, rlat, cufi, cvfi, rugos,         flux_t, flux_q, flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, &
11       .                  debut, lafin, agesno,rugoro,         coefh, t2m, q2m, u10m_srf, v10m_srf, pblh, capcl, oliqcl, cteicl, pblt, &
12       .                  d_t,d_q,d_u,d_v,d_ts,         therm, plcl, fqcalving, ffonte, run_off_lic_0, albsol, sollw, solsw, &
13       .                  flux_t,flux_q,flux_u,flux_v,cdragh,cdragm,         tsol)
14       .                  q2,  
15       .                  dflux_t,dflux_q,      ! From phylmd/clmain.F, version 1.6, 2005/11/16 14:47:19
16       .                  zcoefh,zu1,zv1, t2m, q2m, u10m, v10m,      ! Author: Z. X. Li (LMD/CNRS)
17  cIM cf. AM : pbl      ! Date: Aug. 18th, 1993
18       .                  pblh,capCL,oliqCL,cteiCL,pblT,      ! Objet : interface de couche limite (diffusion verticale)
19       .                  therm,trmb1,trmb2,trmb3,plcl,  
20       .                  fqcalving,ffonte, run_off_lic_0,      ! Tout ce qui a trait aux traceurs est dans "phytrac". Le calcul
21  cIM "slab" ocean      ! de la couche limite pour les traceurs se fait avec "cltrac" et
22       .                  flux_o, flux_g, tslab, seaice)      ! ne tient pas compte de la diff\'erentiation des sous-fractions
23        ! de sol.
24  !  
25  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/clmain.F,v 1.6 2005/11/16 14:47:19 lmdzadmin Exp $      use cdrag_m, only: cdrag
26  !      use clqh_m, only: clqh
27  c      use clvent_m, only: clvent
28  c      use coef_diff_turb_m, only: coef_diff_turb
29  cAA REM:      USE conf_gcm_m, ONLY: lmt_pas
30  cAA-----      USE conf_phys_m, ONLY: iflag_pbl
31  cAA Tout ce qui a trait au traceurs est dans phytrac maintenant      USE dimphy, ONLY: klev, klon
32  cAA pour l'instant le calcul de la couche limite pour les traceurs      USE dimsoil, ONLY: nsoilmx
33  cAA se fait avec cltrac et ne tient pas compte de la differentiation      use hbtm_m, only: hbtm
34  cAA des sous-fraction de sol.      USE histwrite_phy_m, ONLY: histwrite_phy
35  cAA REM bis :      USE indicesol, ONLY: epsfra, is_lic, is_oce, is_sic, is_ter, nbsrf
36  cAA----------      USE interfoce_lim_m, ONLY: interfoce_lim
37  cAA Pour pouvoir extraire les coefficient d'echanges et le vent      use phyetat0_m, only: masque
38  cAA dans la premiere couche, 3 champs supplementaires ont ete crees      use stdlevvar_m, only: stdlevvar
39  cAA zcoefh,zu1 et zv1. Pour l'instant nous avons moyenne les valeurs      USE suphec_m, ONLY: rd, rg, rsigma
40  cAA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir      use time_phylmdz, only: itap
41  cAA si les informations des subsurfaces doivent etre prises en compte  
42  cAA il faudra sortir ces memes champs en leur ajoutant une dimension,      REAL, INTENT(inout):: pctsrf(:, :) ! (klon, nbsrf)
43  cAA c'est a dire nbsrf (nbre de subsurface).      ! pourcentages de surface de chaque maille
44        USE ioipsl  
45        USE interface_surf      REAL, INTENT(IN):: t(klon, klev) ! temperature (K)
46        use dimens_m      REAL, INTENT(IN):: q(klon, klev) ! vapeur d'eau (kg / kg)
47        use indicesol      REAL, INTENT(IN):: u(klon, klev), v(klon, klev) ! vitesse
48        use dimphy      INTEGER, INTENT(IN):: julien ! jour de l'annee en cours
49        use dimsoil      REAL, intent(in):: mu0(klon) ! cosinus de l'angle solaire zenithal    
50        use temps  
51        use iniprint      REAL, INTENT(INout):: ftsol(:, :) ! (klon, nbsrf)
52        use YOMCST      ! skin temperature of surface fraction, in K
53        use yoethf  
54        use fcttre      REAL, INTENT(IN):: cdmmax, cdhmax ! seuils cdrm, cdrh
55        use conf_phys_m  
56        use gath_cpl, only: gath2cpl      REAL, INTENT(inout):: ftsoil(klon, nsoilmx, nbsrf)
57        IMPLICIT none      ! soil temperature of surface fraction
58  c======================================================================  
59  c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818      REAL, INTENT(inout):: qsol(:) ! (klon)
60  c Objet: interface de "couche limite" (diffusion verticale)      ! column-density of water in soil, in kg m-2
61  c Arguments:  
62  c dtime----input-R- interval du temps (secondes)      REAL, INTENT(IN):: paprs(klon, klev + 1) ! pression a intercouche (Pa)
63  c itap-----input-I- numero du pas de temps      REAL, INTENT(IN):: play(klon, klev) ! pression au milieu de couche (Pa)
64  c date0----input-R- jour initial  
65  c t--------input-R- temperature (K)      REAL, INTENT(inout):: fsnow(:, :) ! (klon, nbsrf)
66  c q--------input-R- vapeur d'eau (kg/kg)      ! column-density of mass of snow at the surface, in kg m-2
67  c u--------input-R- vitesse u  
68  c v--------input-R- vitesse v      REAL, INTENT(inout):: fqsurf(klon, nbsrf)
69  c ts-------input-R- temperature du sol (en Kelvin)      REAL, intent(inout):: falbe(klon, nbsrf)
70  c paprs----input-R- pression a intercouche (Pa)  
71  c pplay----input-R- pression au milieu de couche (Pa)      REAL, intent(out):: fluxlat(:, :) ! (klon, nbsrf)
72  c radsol---input-R- flux radiatif net (positif vers le sol) en W/m**2      ! flux de chaleur latente, en W m-2
73  c rlat-----input-R- latitude en degree  
74  c rugos----input-R- longeur de rugosite (en m)      REAL, intent(in):: rain_fall(klon)
75  c cufi-----input-R- resolution des mailles en x (m)      ! liquid water mass flux (kg / m2 / s), positive down
76  c cvfi-----input-R- resolution des mailles en y (m)  
77  c      REAL, intent(in):: snow_fall(klon)
78  c d_t------output-R- le changement pour "t"      ! solid water mass flux (kg / m2 / s), positive down
79  c d_q------output-R- le changement pour "q"  
80  c d_u------output-R- le changement pour "u"      REAL, intent(inout):: frugs(klon, nbsrf) ! longueur de rugosit\'e (en m)
81  c d_v------output-R- le changement pour "v"      real, intent(inout):: agesno(:, :) ! (klon, nbsrf)
82  c d_ts-----output-R- le changement pour "ts"      REAL, INTENT(IN):: rugoro(klon)
83  c flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)  
84  c                    (orientation positive vers le bas)      REAL, intent(out):: d_t(:, :), d_q(:, :) ! (klon, klev)
85  c flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)      ! changement pour t et q
86  c flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal  
87  c flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal      REAL, intent(out):: d_u(klon, klev), d_v(klon, klev)
88  c dflux_t derive du flux sensible      ! changement pour "u" et "v"
89  c dflux_q derive du flux latent  
90  cIM "slab" ocean      REAL, intent(out):: flux_t(klon, nbsrf)
91  c flux_g---output-R-  flux glace (pour OCEAN='slab  ')      ! flux de chaleur sensible (c_p T) (W / m2) (orientation positive
92  c flux_o---output-R-  flux ocean (pour OCEAN='slab  ')      ! vers le bas) à la surface
93  c tslab-in/output-R temperature du slab ocean (en Kelvin) ! uniqmnt pour slab  
94  c seaice---output-R-  glace de mer (kg/m2) (pour OCEAN='slab  ')      REAL, intent(out):: flux_q(klon, nbsrf)
95  ccc      ! flux de vapeur d'eau (kg / m2 / s) à la surface
96  c ffonte----Flux thermique utilise pour fondre la neige  
97  c fqcalving-Flux d'eau "perdue" par la surface et necessaire pour limiter la      REAL, intent(out):: flux_u(:, :), flux_v(:, :) ! (klon, nbsrf)
98  c           hauteur de neige, en kg/m2/s      ! tension du vent (flux turbulent de vent) à la surface, en Pa
99  cAA on rajoute en output yu1 et yv1 qui sont les vents dans  
100  cAA la premiere couche      REAL, INTENT(out):: cdragh(klon), cdragm(klon)
101  cAA ces 4 variables sont maintenant traites dans phytrac      real q2(klon, klev + 1, nbsrf)
102  c itr--------input-I- nombre de traceurs  
103  c tr---------input-R- q. de traceurs      ! Ocean slab:
104  c flux_surf--input-R- flux de traceurs a la surface      REAL, INTENT(out):: dflux_t(klon) ! derive du flux sensible
105  c d_tr-------output-R tendance de traceurs      REAL, INTENT(out):: dflux_q(klon) ! derive du flux latent
106  cIM cf. AM : PBL  
107  c trmb1-------deep_cape      REAL, intent(out):: coefh(:, 2:) ! (klon, 2:klev)
108  c trmb2--------inhibition      ! Pour pouvoir extraire les coefficients d'\'echange, le champ
109  c trmb3-------Point Omega      ! "coefh" a \'et\'e cr\'e\'e. Nous avons moyenn\'e les valeurs de
110  c Cape(klon)-------Cape du thermique      ! ce champ sur les quatre sous-surfaces du mod\`ele.
111  c EauLiq(klon)-------Eau liqu integr du thermique  
112  c ctei(klon)-------Critere d'instab d'entrainmt des nuages de CL      REAL, INTENT(inout):: t2m(klon, nbsrf), q2m(klon, nbsrf)
113  c lcl------- Niveau de condensation  
114  c pblh------- HCL      REAL, INTENT(inout):: u10m_srf(:, :), v10m_srf(:, :) ! (klon, nbsrf)
115  c pblT------- T au nveau HCL      ! composantes du vent \`a 10m sans spirale d'Ekman
116  c======================================================================  
117  c$$$ PB ajout pour soil      ! Ionela Musat. Cf. Anne Mathieu : planetary boundary layer, hbtm.
118  c      ! Comme les autres diagnostics on cumule dans physiq ce qui permet
119        REAL dtime      ! de sortir les grandeurs par sous-surface.
120        real date0      REAL pblh(klon, nbsrf) ! height of planetary boundary layer
121        integer itap      REAL capcl(klon, nbsrf)
122        REAL t(klon,klev), q(klon,klev)      REAL oliqcl(klon, nbsrf)
123        REAL u(klon,klev), v(klon,klev)      REAL cteicl(klon, nbsrf)
124  cIM 230604 BAD  REAL radsol(klon) ???      REAL, INTENT(inout):: pblt(klon, nbsrf) ! T au nveau HCL
125        REAL, intent(in):: paprs(klon,klev+1)      REAL therm(klon, nbsrf)
126        real pplay(klon,klev)      REAL plcl(klon, nbsrf)
127        REAL, intent(in):: rlon(klon), rlat(klon)  
128        real cufi(klon), cvfi(klon)      REAL, intent(out):: fqcalving(klon, nbsrf)
129        REAL d_t(klon, klev), d_q(klon, klev)      ! flux d'eau "perdue" par la surface et necessaire pour limiter la
130        REAL d_u(klon, klev), d_v(klon, klev)      ! hauteur de neige, en kg / m2 / s
131        REAL flux_t(klon,klev, nbsrf), flux_q(klon,klev, nbsrf)  
132        REAL dflux_t(klon), dflux_q(klon)      real ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
133  cIM "slab" ocean      REAL, intent(inout):: run_off_lic_0(:) ! (klon)
134        REAL flux_o(klon), flux_g(klon)  
135        REAL y_flux_o(klon), y_flux_g(klon)      REAL, intent(out):: albsol(:) ! (klon)
136        REAL tslab(klon), ytslab(klon)      ! albedo du sol total, visible, moyen par maille
137        REAL seaice(klon), y_seaice(klon)  
138  cIM cf JLD      REAL, intent(in):: sollw(:) ! (klon)
139        REAL y_fqcalving(klon), y_ffonte(klon)      ! surface net downward longwave flux, in W m-2
140        REAL fqcalving(klon,nbsrf), ffonte(klon,nbsrf)  
141        REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      REAL, intent(in):: solsw(:) ! (klon)
142        ! surface net downward shortwave flux, in W m-2
143        REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)  
144        REAL rugmer(klon), agesno(klon,nbsrf),rugoro(klon)      REAL, intent(in):: tsol(:) ! (klon)
145        REAL cdragh(klon), cdragm(klon)  
146        integer jour            ! jour de l'annee en cours      ! Local:
147        real rmu0(klon)         ! cosinus de l'angle solaire zenithal  
148        REAL co2_ppm            ! taux CO2 atmosphere      REAL d_ts(klon, nbsrf) ! variation of ftsol
149        LOGICAL, intent(in):: debut      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
150        logical, intent(in):: lafin      REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
151        logical ok_veget  
152        character*6 ocean      ! la nouvelle repartition des surfaces sortie de l'interface
153        integer npas, nexca      REAL, save:: pctsrf_new_oce(klon)
154  c      REAL, save:: pctsrf_new_sic(klon)
155        REAL pctsrf(klon,nbsrf)  
156        REAL ts(klon,nbsrf)      REAL y_fqcalving(klon), y_ffonte(klon)
157        REAL d_ts(klon,nbsrf)      real y_run_off_lic_0(klon), y_run_off_lic(klon)
158        REAL snow(klon,nbsrf)      REAL run_off_lic(klon) ! ruissellement total
159        REAL qsurf(klon,nbsrf)      REAL rugmer(klon)
160        REAL evap(klon,nbsrf)      REAL ytsoil(klon, nsoilmx)
161        REAL albe(klon,nbsrf)      REAL yts(klon), ypctsrf(klon), yz0_new(klon)
162        REAL alblw(klon,nbsrf)      real yrugos(klon) ! longueur de rugosite (en m)
163  c$$$ PB      REAL yalb(klon)
164        REAL fluxlat(klon,nbsrf)      REAL snow(klon) ! column-density of mass of snow at the surface, in kg m-2
165  C      real yqsurf(klon), yagesno(klon)
166        real rain_f(klon), snow_f(klon)      real yqsol(klon) ! column-density of water in soil, in kg m-2
167        REAL fder(klon)      REAL yrain_fall(klon) ! liquid water mass flux (kg / m2 / s), positive down
168  cIM cf. JLD   REAL sollw(klon), solsw(klon), sollwdown(klon)      REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
169        REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)      REAL yrugm(klon), radsol(klon), yrugoro(klon)
170        REAL rugos(klon,nbsrf)      REAL yfluxlat(klon)
171  C la nouvelle repartition des surfaces sortie de l'interface      REAL y_d_ts(klon)
172        REAL pctsrf_new(klon,nbsrf)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
173  cAA      REAL y_d_u(klon, klev), y_d_v(klon, klev)
174        REAL zcoefh(klon,klev)      REAL y_flux_t(klon), y_flux_q(klon)
175        REAL zu1(klon)      REAL y_flux_u(klon), y_flux_v(klon)
176        REAL zv1(klon)      REAL y_dflux_t(klon), y_dflux_q(klon)
177  cAA      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
178  c$$$ PB ajout pour soil      real ycdragh(klon), ycdragm(klon)
179        LOGICAL soil_model      REAL yu(klon, klev), yv(klon, klev)
180  cIM ajout seuils cdrm, cdrh      REAL yt(klon, klev), yq(klon, klev)
181        REAL cdmmax, cdhmax      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
182  cIM: 261103      REAL yq2(klon, klev + 1)
183        REAL ksta, ksta_ter      REAL delp(klon, klev)
184        LOGICAL ok_kzmin      INTEGER i, k, nsrf
185  cIM: 261103      INTEGER ni(klon), knon, j
186        REAL ftsoil(klon,nsoilmx,nbsrf)  
187        REAL ytsoil(klon,nsoilmx)      REAL pctsrf_pot(klon, nbsrf)
188        REAL qsol(klon)      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
189  c======================================================================      ! apparitions ou disparitions de la glace de mer
190        EXTERNAL clqh, clvent, coefkz, calbeta, cltrac  
191  c======================================================================      REAL yt2m(klon), yq2m(klon), wind10m(klon)
192        REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL ustar(klon)
193        REAL yalb(klon)  
194        REAL yalblw(klon)      REAL yt10m(klon), yq10m(klon)
195        REAL yu1(klon), yv1(klon)      REAL ypblh(klon)
196        real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL ylcl(klon)
197        real yrain_f(klon), ysnow_f(klon)      REAL ycapcl(klon)
198        real ysollw(klon), ysolsw(klon), ysollwdown(klon)      REAL yoliqcl(klon)
199        real yfder(klon), ytaux(klon), ytauy(klon)      REAL ycteicl(klon)
200        REAL yrugm(klon), yrads(klon),yrugoro(klon)      REAL ypblt(klon)
201  c$$$ PB      REAL ytherm(klon)
202        REAL yfluxlat(klon)      REAL u1(klon), v1(klon)
203  C      REAL tair1(klon), qair1(klon), tairsol(klon)
204        REAL y_d_ts(klon)      REAL psfce(klon), patm(klon)
205        REAL y_d_t(klon, klev), y_d_q(klon, klev)      REAL zgeo1(klon)
206        REAL y_d_u(klon, klev), y_d_v(klon, klev)      REAL rugo1(klon)
207        REAL y_flux_t(klon,klev), y_flux_q(klon,klev)      REAL zgeop(klon, klev)
208        REAL y_flux_u(klon,klev), y_flux_v(klon,klev)  
209        REAL y_dflux_t(klon), y_dflux_q(klon)      !------------------------------------------------------------
210        REAL ycoefh(klon,klev), ycoefm(klon,klev)  
211        REAL yu(klon,klev), yv(klon,klev)      albsol = sum(falbe * pctsrf, dim = 2)
212        REAL yt(klon,klev), yq(klon,klev)  
213        REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)      ! R\'epartition sous maille des flux longwave et shortwave
214  c      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
215        LOGICAL ok_nonloc  
216        PARAMETER (ok_nonloc=.FALSE.)      forall (nsrf = 1:nbsrf)
217        REAL ycoefm0(klon,klev), ycoefh0(klon,klev)         fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
218                * (tsol - ftsol(:, nsrf))
219  cIM 081204 hcl_Anne ? BEG         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
220        real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)      END forall
221        real ykmm(klon,klev+1),ykmn(klon,klev+1)  
222        real ykmq(klon,klev+1)      ytherm = 0.
223        real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)  
224        real q2diag(klon,klev+1)      DO k = 1, klev ! epaisseur de couche
225  cIM 081204   real yustar(klon),y_cd_m(klon),y_cd_h(klon)         DO i = 1, klon
226  cIM 081204 hcl_Anne ? END            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
227  c         END DO
228        REAL u1lay(klon), v1lay(klon)      END DO
229        REAL delp(klon,klev)  
230        INTEGER i, k, nsrf      ! Initialization:
231  cAA   INTEGER it      rugmer = 0.
232        INTEGER ni(klon), knon, j      cdragh = 0.
233  c Introduction d'une variable "pourcentage potentiel" pour tenir compte      cdragm = 0.
234  c des eventuelles apparitions et/ou disparitions de la glace de mer      dflux_t = 0.
235        REAL pctsrf_pot(klon,nbsrf)      dflux_q = 0.
236              yrugos = 0.
237  c======================================================================      ypaprs = 0.
238        REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      ypplay = 0.
239  c======================================================================      ydelp = 0.
240  c      yrugoro = 0.
241  c maf pour sorties IOISPL en cas de debugagage      d_ts = 0.
242  c      flux_t = 0.
243        CHARACTER*80 cldebug      flux_q = 0.
244        SAVE cldebug      flux_u = 0.
245        CHARACTER*8 cl_surf(nbsrf)      flux_v = 0.
246        SAVE cl_surf      fluxlat = 0.
247        INTEGER nhoridbg, nidbg      d_t = 0.
248        SAVE nhoridbg, nidbg      d_q = 0.
249        INTEGER ndexbg(iim*(jjm+1))      d_u = 0.
250        REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian      d_v = 0.
251        REAL tabindx(klon)      coefh = 0.
252        REAL debugtab(iim,jjm+1)      fqcalving = 0.
253        LOGICAL first_appel      run_off_lic = 0.
254        SAVE first_appel  
255        DATA first_appel/.true./      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
256        LOGICAL debugindex      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
257        SAVE debugindex      ! (\`a affiner).
258        DATA debugindex/.false./  
259        integer idayref      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
260        REAL t2m(klon,nbsrf), q2m(klon,nbsrf)      pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
261        REAL u10m(klon,nbsrf), v10m(klon,nbsrf)      pctsrf_pot(:, is_oce) = 1. - masque
262  c      pctsrf_pot(:, is_sic) = 1. - masque
263        REAL yt2m(klon), yq2m(klon), yu10m(klon)  
264        REAL yustar(klon)      ! Tester si c'est le moment de lire le fichier:
265  c -- LOOP      if (mod(itap - 1, lmt_pas) == 0) then
266         REAL yu10mx(klon)         CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
267         REAL yu10my(klon)      endif
268         REAL ywindsp(klon)  
269  c -- LOOP      ! Boucler sur toutes les sous-fractions du sol:
270  c  
271        REAL yt10m(klon), yq10m(klon)      loop_surface: DO nsrf = 1, nbsrf
272  cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui         ! Define ni and knon:
273  c   permet de sortir les grdeurs par sous surface)  
274        REAL pblh(klon,nbsrf)         ni = 0
275        REAL plcl(klon,nbsrf)         knon = 0
276        REAL capCL(klon,nbsrf)  
277        REAL oliqCL(klon,nbsrf)         DO i = 1, klon
278        REAL cteiCL(klon,nbsrf)            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
279        REAL pblT(klon,nbsrf)            ! "potentielles"
280        REAL therm(klon,nbsrf)            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
281        REAL trmb1(klon,nbsrf)               knon = knon + 1
282        REAL trmb2(klon,nbsrf)               ni(knon) = i
283        REAL trmb3(klon,nbsrf)            END IF
284        REAL ypblh(klon)         END DO
285        REAL ylcl(klon)  
286        REAL ycapCL(klon)         if_knon: IF (knon /= 0) then
287        REAL yoliqCL(klon)            ypctsrf(:knon) = pctsrf(ni(:knon), nsrf)
288        REAL ycteiCL(klon)            yts(:knon) = ftsol(ni(:knon), nsrf)
289        REAL ypblT(klon)            snow(:knon) = fsnow(ni(:knon), nsrf)
290        REAL ytherm(klon)            yqsurf(:knon) = fqsurf(ni(:knon), nsrf)
291        REAL ytrmb1(klon)            yalb(:knon) = falbe(ni(:knon), nsrf)
292        REAL ytrmb2(klon)            yrain_fall(:knon) = rain_fall(ni(:knon))
293        REAL ytrmb3(klon)            ysnow_fall(:knon) = snow_fall(ni(:knon))
294        REAL y_cd_h(klon), y_cd_m(klon)            yagesno(:knon) = agesno(ni(:knon), nsrf)
295  c     REAL ygamt(klon,2:klev) ! contre-gradient pour temperature            yrugos(:knon) = frugs(ni(:knon), nsrf)
296  c     REAL ygamq(klon,2:klev) ! contre-gradient pour humidite            yrugoro(:knon) = rugoro(ni(:knon))
297        REAL uzon(klon), vmer(klon)            radsol(:knon) = fsolsw(ni(:knon), nsrf) + fsollw(ni(:knon), nsrf)
298        REAL tair1(klon), qair1(klon), tairsol(klon)            ypaprs(:knon, klev + 1) = paprs(ni(:knon), klev + 1)
299        REAL psfce(klon), patm(klon)            y_run_off_lic_0(:knon) = run_off_lic_0(ni(:knon))
300  c  
301        REAL qairsol(klon), zgeo1(klon)            ! For continent, copy soil water content
302        REAL rugo1(klon)            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
303  c  
304        LOGICAL zxli ! utiliser un jeu de fonctions simples            ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
305        PARAMETER (zxli=.FALSE.)  
306  c            DO k = 1, klev
307        REAL zt, zqs, zdelta, zcor               DO j = 1, knon
308        REAL t_coup                  i = ni(j)
309        PARAMETER(t_coup=273.15)                  ypaprs(j, k) = paprs(i, k)
310  C                  ypplay(j, k) = play(i, k)
311        character (len = 20) :: modname = 'clmain'                  ydelp(j, k) = delp(i, k)
312        LOGICAL check                  yu(j, k) = u(i, k)
313        PARAMETER (check=.false.)                  yv(j, k) = v(i, k)
314                    yt(j, k) = t(i, k)
315                    yq(j, k) = q(i, k)
316  c initialisation Anne               END DO
317        ytherm(:) = 0.            END DO
318  C  
319        if (check) THEN            ! Calculer les géopotentiels de chaque couche:
320            write(*,*) modname,'  klon=',klon  
321  CC        call flush(6)            zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
322        endif                 + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
323        IF (debugindex .and. first_appel) THEN  
324            first_appel=.false.            DO k = 2, klev
325  !               zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
326  ! initialisation sorties netcdf                    * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
327  !                    * (ypplay(:knon, k - 1) - ypplay(:knon, k))
           idayref = day_ini  
           CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)  
           CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)  
           DO i = 1, iim  
             zx_lon(i,1) = rlon(i+1)  
             zx_lon(i,jjm+1) = rlon(i+1)  
328            ENDDO            ENDDO
329            CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)  
330            cldebug='sous_index'            CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
331            CALL histbeg_totreg(cldebug, iim,zx_lon(:,1),jjm+1,                 yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
332       $        zx_lat(1,:),1,iim,1,jjm                 yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
333       $        +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)                 ycdragh(:knon))
334  ! no vertical axis  
335            cl_surf(1)='ter'            IF (iflag_pbl == 1) THEN
336            cl_surf(2)='lic'               ycdragm(:knon) = max(ycdragm(:knon), 0.)
337            cl_surf(3)='oce'               ycdragh(:knon) = max(ycdragh(:knon), 0.)
338            cl_surf(4)='sic'            end IF
339            DO nsrf=1,nbsrf  
340              CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,            ! on met un seuil pour ycdragm et ycdragh
341       $          jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)            IF (nsrf == is_oce) THEN
342                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
343                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
344              END IF
345    
346              IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
347              call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
348                   ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
349                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
350                   ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
351    
352              CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
353                   ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
354                   ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
355                   y_flux_u(:knon))
356              CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
357                   ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
358                   ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
359                   y_flux_v(:knon))
360    
361              CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
362                   mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
363                   yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
364                   yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
365                   ydelp(:knon, :), radsol(:knon), yalb(:knon), snow(:knon), &
366                   yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
367                   yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
368                   y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
369                   yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
370                   y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
371                   y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
372    
373              ! calculer la longueur de rugosite sur ocean
374    
375              yrugm = 0.
376    
377              IF (nsrf == is_oce) THEN
378                 DO j = 1, knon
379                    yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
380                         / rg + 0.11 * 14E-6 &
381                         / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
382                    yrugm(j) = max(1.5E-05, yrugm(j))
383                 END DO
384              END IF
385    
386              DO k = 1, klev
387                 DO j = 1, knon
388                    i = ni(j)
389                    y_d_t(j, k) = y_d_t(j, k) * ypctsrf(j)
390                    y_d_q(j, k) = y_d_q(j, k) * ypctsrf(j)
391                    y_d_u(j, k) = y_d_u(j, k) * ypctsrf(j)
392                    y_d_v(j, k) = y_d_v(j, k) * ypctsrf(j)
393                 END DO
394              END DO
395    
396              flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
397              flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
398              flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
399              flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
400    
401              falbe(:, nsrf) = 0.
402              fsnow(:, nsrf) = 0.
403              fqsurf(:, nsrf) = 0.
404              frugs(:, nsrf) = 0.
405              DO j = 1, knon
406                 i = ni(j)
407                 d_ts(i, nsrf) = y_d_ts(j)
408                 falbe(i, nsrf) = yalb(j)
409                 fsnow(i, nsrf) = snow(j)
410                 fqsurf(i, nsrf) = yqsurf(j)
411                 frugs(i, nsrf) = yz0_new(j)
412                 fluxlat(i, nsrf) = yfluxlat(j)
413                 IF (nsrf == is_oce) THEN
414                    rugmer(i) = yrugm(j)
415                    frugs(i, nsrf) = yrugm(j)
416                 END IF
417                 agesno(i, nsrf) = yagesno(j)
418                 fqcalving(i, nsrf) = y_fqcalving(j)
419                 ffonte(i, nsrf) = y_ffonte(j)
420                 cdragh(i) = cdragh(i) + ycdragh(j) * ypctsrf(j)
421                 cdragm(i) = cdragm(i) + ycdragm(j) * ypctsrf(j)
422                 dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypctsrf(j)
423                 dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypctsrf(j)
424              END DO
425              IF (nsrf == is_ter) THEN
426                 qsol(ni(:knon)) = yqsol(:knon)
427              else IF (nsrf == is_lic) THEN
428                 DO j = 1, knon
429                    i = ni(j)
430                    run_off_lic_0(i) = y_run_off_lic_0(j)
431                    run_off_lic(i) = y_run_off_lic(j)
432                 END DO
433              END IF
434    
435              ftsoil(:, :, nsrf) = 0.
436              ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
437    
438              DO j = 1, knon
439                 i = ni(j)
440                 DO k = 1, klev
441                    d_t(i, k) = d_t(i, k) + y_d_t(j, k)
442                    d_q(i, k) = d_q(i, k) + y_d_q(j, k)
443                    d_u(i, k) = d_u(i, k) + y_d_u(j, k)
444                    d_v(i, k) = d_v(i, k) + y_d_v(j, k)
445                 END DO
446              END DO
447    
448              forall (k = 2:klev) coefh(ni(:knon), k) &
449                   = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypctsrf(:knon)
450    
451              ! diagnostic t, q a 2m et u, v a 10m
452    
453              DO j = 1, knon
454                 i = ni(j)
455                 u1(j) = yu(j, 1) + y_d_u(j, 1)
456                 v1(j) = yv(j, 1) + y_d_v(j, 1)
457                 tair1(j) = yt(j, 1) + y_d_t(j, 1)
458                 qair1(j) = yq(j, 1) + y_d_q(j, 1)
459                 zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
460                      1))) * (ypaprs(j, 1)-ypplay(j, 1))
461                 tairsol(j) = yts(j) + y_d_ts(j)
462                 rugo1(j) = yrugos(j)
463                 IF (nsrf == is_oce) THEN
464                    rugo1(j) = frugs(i, nsrf)
465                 END IF
466                 psfce(j) = ypaprs(j, 1)
467                 patm(j) = ypplay(j, 1)
468              END DO
469    
470              CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
471                   zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
472                   yt10m, yq10m, wind10m(:knon), ustar(:knon))
473    
474              DO j = 1, knon
475                 i = ni(j)
476                 t2m(i, nsrf) = yt2m(j)
477                 q2m(i, nsrf) = yq2m(j)
478    
479                 u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
480                      / sqrt(u1(j)**2 + v1(j)**2)
481                 v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
482                      / sqrt(u1(j)**2 + v1(j)**2)
483            END DO            END DO
484            CALL histend(nidbg)  
485            CALL histsync(nidbg)            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
486        ENDIF                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
487                             yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
488        DO k = 1, klev   ! epaisseur de couche                 ytherm, ylcl)
489        DO i = 1, klon  
          delp(i,k) = paprs(i,k)-paprs(i,k+1)  
       ENDDO  
       ENDDO  
       DO i = 1, klon  ! vent de la premiere couche  
 ccc         zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))  
          zx_alf1 = 1.0  
          zx_alf2 = 1.0 - zx_alf1  
          u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2  
          v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2  
       ENDDO  
 c  
 c initialisation:  
 c  
       DO i = 1, klon  
          rugmer(i) = 0.0  
          cdragh(i) = 0.0  
          cdragm(i) = 0.0  
          dflux_t(i) = 0.0  
          dflux_q(i) = 0.0  
          zu1(i) = 0.0  
          zv1(i) = 0.0  
       ENDDO  
       ypct = 0.0  
       yts = 0.0  
       ysnow = 0.0  
       yqsurf = 0.0  
       yalb = 0.0  
       yalblw = 0.0  
       yrain_f = 0.0  
       ysnow_f = 0.0  
       yfder = 0.0  
       ytaux = 0.0  
       ytauy = 0.0  
       ysolsw = 0.0  
       ysollw = 0.0  
       ysollwdown = 0.0  
       yrugos = 0.0  
       yu1 = 0.0  
       yv1 = 0.0  
       yrads = 0.0  
       ypaprs = 0.0  
       ypplay = 0.0  
       ydelp = 0.0  
       yu = 0.0  
       yv = 0.0  
       yt = 0.0  
       yq = 0.0  
       pctsrf_new = 0.0  
       y_flux_u = 0.0  
       y_flux_v = 0.0  
 C$$ PB  
       y_dflux_t = 0.0  
       y_dflux_q = 0.0  
       ytsoil = 999999.  
       yrugoro = 0.  
 c -- LOOP  
       yu10mx = 0.0  
       yu10my = 0.0  
       ywindsp = 0.0  
 c -- LOOP  
       DO nsrf = 1, nbsrf  
       DO i = 1, klon  
          d_ts(i,nsrf) = 0.0  
       ENDDO  
       END DO  
 C§§§ PB  
       yfluxlat=0.  
       flux_t = 0.  
       flux_q = 0.  
       flux_u = 0.  
       flux_v = 0.  
       DO k = 1, klev  
       DO i = 1, klon  
          d_t(i,k) = 0.0  
          d_q(i,k) = 0.0  
 c$$$         flux_t(i,k) = 0.0  
 c$$$         flux_q(i,k) = 0.0  
          d_u(i,k) = 0.0  
          d_v(i,k) = 0.0  
 c$$$         flux_u(i,k) = 0.0  
 c$$$         flux_v(i,k) = 0.0  
          zcoefh(i,k) = 0.0  
       ENDDO  
       ENDDO  
 cAA      IF (itr.GE.1) THEN  
 cAA      DO it = 1, itr  
 cAA      DO k = 1, klev  
 cAA      DO i = 1, klon  
 cAA         d_tr(i,k,it) = 0.0  
 cAA      ENDDO  
 cAA      ENDDO  
 cAA      ENDDO  
 cAA      ENDIF  
   
 c  
 c Boucler sur toutes les sous-fractions du sol:  
 c  
 C Initialisation des "pourcentages potentiels". On considere ici qu'on  
 C peut avoir potentiellementdela glace sur tout le domaine oceanique  
 C (a affiner)  
   
       pctsrf_pot = pctsrf  
       pctsrf_pot(:,is_oce) = 1. - zmasq(:)  
       pctsrf_pot(:,is_sic) = 1. - zmasq(:)  
   
       DO 99999 nsrf = 1, nbsrf  
   
 c chercher les indices:  
       DO j = 1, klon  
          ni(j) = 0  
       ENDDO  
       knon = 0  
       DO i = 1, klon  
   
 C pour determiner le domaine a traiter on utilise les surfaces "potentielles"  
 C    
       IF (pctsrf_pot(i,nsrf).GT.epsfra) THEN  
          knon = knon + 1  
          ni(knon) = i  
       ENDIF  
       ENDDO  
 c  
       if (check) THEN  
           write(*,*)'CLMAIN, nsrf, knon =',nsrf, knon  
 CC        call flush(6)  
       endif  
 c  
 c variables pour avoir une sortie IOIPSL des INDEX  
 c  
       IF (debugindex) THEN  
           tabindx(:)=0.  
 c          tabindx(1:knon)=(/FLOAT(i),i=1:knon/)  
           DO i=1,knon  
             tabindx(1:knon)=FLOAT(i)  
           END DO  
           debugtab(:,:)=0.  
           ndexbg(:)=0  
           CALL gath2cpl(tabindx,debugtab,klon,knon,iim,jjm,ni)  
           CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1)  
      $        ,ndexbg)  
       ENDIF  
       IF (knon.EQ.0) GOTO 99999  
       DO j = 1, knon  
       i = ni(j)  
         ypct(j) = pctsrf(i,nsrf)  
         yts(j) = ts(i,nsrf)  
 cIM "slab" ocean  
 c        PRINT *, 'tslab = ', i, tslab(i)  
         ytslab(i) = tslab(i)  
 c  
         ysnow(j) = snow(i,nsrf)  
         yqsurf(j) = qsurf(i,nsrf)  
         yalb(j) = albe(i,nsrf)  
         yalblw(j) = alblw(i,nsrf)  
         yrain_f(j) = rain_f(i)  
         ysnow_f(j) = snow_f(i)  
         yagesno(j) = agesno(i,nsrf)  
         yfder(j) = fder(i)  
         ytaux(j) = flux_u(i,1,nsrf)  
         ytauy(j) = flux_v(i,1,nsrf)  
         ysolsw(j) = solsw(i,nsrf)  
         ysollw(j) = sollw(i,nsrf)  
         ysollwdown(j) = sollwdown(i)  
         yrugos(j) = rugos(i,nsrf)  
         yrugoro(j) = rugoro(i)  
         yu1(j) = u1lay(i)  
         yv1(j) = v1lay(i)  
         yrads(j) =  ysolsw(j)+ ysollw(j)  
         ypaprs(j,klev+1) = paprs(i,klev+1)  
         y_run_off_lic_0(j) = run_off_lic_0(i)  
 c -- LOOP  
        yu10mx(j) = u10m(i,nsrf)  
        yu10my(j) = v10m(i,nsrf)  
        ywindsp(j) = SQRT(yu10mx(j)*yu10mx(j) + yu10my(j)*yu10my(j) )  
 c -- LOOP  
       END DO  
 C  
 C     IF bucket model for continent, copy soil water content  
       IF ( nsrf .eq. is_ter .and. .not. ok_veget ) THEN  
490            DO j = 1, knon            DO j = 1, knon
491              i = ni(j)               i = ni(j)
492              yqsol(j) = qsol(i)               pblh(i, nsrf) = ypblh(j)
493                 plcl(i, nsrf) = ylcl(j)
494                 capcl(i, nsrf) = ycapcl(j)
495                 oliqcl(i, nsrf) = yoliqcl(j)
496                 cteicl(i, nsrf) = ycteicl(j)
497                 pblt(i, nsrf) = ypblt(j)
498                 therm(i, nsrf) = ytherm(j)
499            END DO            END DO
       ELSE  
           yqsol(:)=0.  
       ENDIF  
 c$$$ PB ajour pour soil  
       DO k = 1, nsoilmx  
         DO j = 1, knon  
           i = ni(j)  
           ytsoil(j,k) = ftsoil(i,k,nsrf)  
         END DO    
       END DO  
       DO k = 1, klev  
       DO j = 1, knon  
       i = ni(j)  
         ypaprs(j,k) = paprs(i,k)  
         ypplay(j,k) = pplay(i,k)  
         ydelp(j,k) = delp(i,k)  
         yu(j,k) = u(i,k)  
         yv(j,k) = v(i,k)  
         yt(j,k) = t(i,k)  
         yq(j,k) = q(i,k)  
       ENDDO  
       ENDDO  
 c  
 c  
 c calculer Cdrag et les coefficients d'echange  
       CALL coefkz(nsrf, knon, ypaprs, ypplay,  
 cIM 261103  
      .     ksta, ksta_ter,  
 cIM 261103  
      .            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)  
         DO k = 1, klev  
         DO i = 1, knon  
            ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))  
            ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))  
         ENDDO  
         ENDDO  
       endif  
 c  
 cIM cf JLD : on seuille ycoefm et ycoefh  
       if (nsrf.eq.is_oce) then  
          do j=1,knon  
 c           ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)  
             ycoefm(j,1)=min(ycoefm(j,1),cdmmax)  
 c           ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)  
             ycoefh(j,1)=min(ycoefh(j,1),cdhmax)  
          enddo  
       endif  
   
 c  
 cIM: 261103  
       if (ok_kzmin) THEN  
 cIM cf FH: 201103 BEG  
 c   Calcul d'une diffusion minimale pour les conditions tres stables.  
       call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm  
      .   ,ycoefm0,ycoefh0)  
 c      call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ         ')  
 c      call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN      ')  
   
        if ( 1.eq.1 ) then  
        DO k = 1, klev  
        DO i = 1, knon  
           ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))  
           ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))  
        ENDDO  
        ENDDO  
        endif  
 cIM cf FH: 201103 END  
       endif !ok_kzmin  
 cIM: 261103  
   
   
       IF (iflag_pbl.ge.3) then  
   
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
 c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin  
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
   
          yzlay(1:knon,1)=  
      .   RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1)))  
      .   *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG  
          do k=2,klev  
             yzlay(1:knon,k)=  
      .      yzlay(1:knon,k-1)+RD*0.5*(yt(1:knon,k-1)+yt(1:knon,k))  
      .      /ypaprs(1:knon,k)*(ypplay(1:knon,k-1)-ypplay(1:knon,k))/RG  
          enddo  
          do k=1,klev  
             yteta(1:knon,k)=  
      .      yt(1:knon,k)*(ypaprs(1:knon,1)/ypplay(1:knon,k))**rkappa  
      .      *(1.+0.61*yq(1:knon,k))  
          enddo  
          yzlev(1:knon,1)=0.  
          yzlev(1:knon,klev+1)=2.*yzlay(1:knon,klev)-yzlay(1:knon,klev-1)  
          do k=2,klev  
             yzlev(1:knon,k)=0.5*(yzlay(1:knon,k)+yzlay(1:knon,k-1))  
          enddo  
          DO k = 1, klev+1  
             DO j = 1, knon  
                i = ni(j)  
                yq2(j,k)=q2(i,k,nsrf)  
             enddo  
          enddo  
   
   
 c   Bug introduit volontairement pour converger avec les resultats  
 c  du papier sur les thermiques.  
          if (1.eq.1) then  
          y_cd_m(1:knon) = ycoefm(1:knon,1)  
          y_cd_h(1:knon) = ycoefh(1:knon,1)  
          else  
          y_cd_h(1:knon) = ycoefm(1:knon,1)  
          y_cd_m(1:knon) = ycoefh(1:knon,1)  
          endif  
          call ustarhb(knon,yu,yv,y_cd_m, yustar)  
   
         if (prt_level > 9) THEN  
           WRITE(lunout,*)'USTAR = ',yustar  
         ENDIF  
   
 c   iflag_pbl peut etre utilise comme longuer de melange  
   
          if (iflag_pbl.ge.11) then  
             call vdif_kcay(knon,dtime,rg,rd,ypaprs,yt  
      s      ,yzlev,yzlay,yu,yv,yteta  
      s      ,y_cd_m,yq2,q2diag,ykmm,ykmn,yustar,  
      s      iflag_pbl)  
          else  
             call yamada4(knon,dtime,rg,rd,ypaprs,yt  
      s      ,yzlev,yzlay,yu,yv,yteta  
      s      ,y_cd_m,yq2,ykmm,ykmn,ykmq,yustar,  
      s      iflag_pbl)  
          endif  
   
          ycoefm(1:knon,1)=y_cd_m(1:knon)  
          ycoefh(1:knon,1)=y_cd_h(1:knon)  
          ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)  
          ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)  
   
   
       ENDIF  
   
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
 c calculer la diffusion des vitesses "u" et "v"  
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
   
       CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,  
      s            y_d_u,y_flux_u)  
       CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,  
      s            y_d_v,y_flux_v)  
   
 c pour le couplage  
       ytaux = y_flux_u(:,1)  
       ytauy = y_flux_v(:,1)  
   
 c FH modif sur le cdrag temperature  
 c$$$PB : déplace dans clcdrag  
 c$$$      do i=1,knon  
 c$$$         ycoefh(i,1)=ycoefm(i,1)*0.8  
 c$$$      enddo  
   
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
 c calculer la diffusion de "q" et de "h"  
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
       CALL clqh(dtime, itap, date0,jour, debut,lafin,  
      e          rlon, rlat, cufi, cvfi,  
      e          knon, nsrf, ni, pctsrf,  
      e          soil_model, ytsoil,yqsol,  
      e          ok_veget, ocean, npas, nexca,  
      e          rmu0, co2_ppm, yrugos, yrugoro,  
      e          yu1, yv1, ycoefh,  
      e          yt,yq,yts,ypaprs,ypplay,  
      e          ydelp,yrads,yalb, yalblw, ysnow, yqsurf,  
      e          yrain_f, ysnow_f, yfder, ytaux, ytauy,  
 c -- LOOP  
      e          ywindsp,  
 c -- LOOP  
 c$$$     e          ysollw, ysolsw,  
      e          ysollw, ysollwdown, ysolsw,yfluxlat,  
      s          pctsrf_new, yagesno,  
      s          y_d_t, y_d_q, y_d_ts, yz0_new,  
      s          y_flux_t, y_flux_q, y_dflux_t, y_dflux_q,  
      s          y_fqcalving,y_ffonte,y_run_off_lic_0,  
 cIM "slab" ocean  
      s          y_flux_o, y_flux_g, ytslab, y_seaice)  
 c  
 c calculer la longueur de rugosite sur ocean  
       yrugm=0.  
       IF (nsrf.EQ.is_oce) THEN  
       DO j = 1, knon  
          yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG  
      $      +  0.11*14e-6 / sqrt(ycoefm(j,1) * (yu1(j)**2+yv1(j)**2))  
          yrugm(j) = MAX(1.5e-05,yrugm(j))  
       ENDDO  
       ENDIF  
       DO j = 1, knon  
          y_dflux_t(j) = y_dflux_t(j) * ypct(j)  
          y_dflux_q(j) = y_dflux_q(j) * ypct(j)  
          yu1(j) = yu1(j) *  ypct(j)  
          yv1(j) = yv1(j) *  ypct(j)  
       ENDDO  
 c  
       DO k = 1, klev  
         DO j = 1, knon  
           i = ni(j)  
           ycoefh(j,k) = ycoefh(j,k) * ypct(j)  
           ycoefm(j,k) = ycoefm(j,k) * ypct(j)  
           y_d_t(j,k) = y_d_t(j,k) * ypct(j)  
           y_d_q(j,k) = y_d_q(j,k) * ypct(j)  
 C§§§ PB  
           flux_t(i,k,nsrf) = y_flux_t(j,k)  
           flux_q(i,k,nsrf) = y_flux_q(j,k)  
           flux_u(i,k,nsrf) = y_flux_u(j,k)  
           flux_v(i,k,nsrf) = y_flux_v(j,k)  
 c$$$ PB        y_flux_t(j,k) = y_flux_t(j,k) * ypct(j)  
 c$$$ PB        y_flux_q(j,k) = y_flux_q(j,k) * ypct(j)  
           y_d_u(j,k) = y_d_u(j,k) * ypct(j)  
           y_d_v(j,k) = y_d_v(j,k) * ypct(j)  
 c$$$ PB        y_flux_u(j,k) = y_flux_u(j,k) * ypct(j)  
 c$$$ PB        y_flux_v(j,k) = y_flux_v(j,k) * ypct(j)  
         ENDDO  
       ENDDO  
   
   
       evap(:,nsrf) = - flux_q(:,1,nsrf)  
 c  
       albe(:, nsrf) = 0.  
       alblw(:, nsrf) = 0.  
       snow(:, nsrf) = 0.  
       qsurf(:, nsrf) = 0.  
       rugos(:, nsrf) = 0.  
       fluxlat(:,nsrf) = 0.  
       DO j = 1, knon  
          i = ni(j)  
          d_ts(i,nsrf) = y_d_ts(j)  
          albe(i,nsrf) = yalb(j)  
          alblw(i,nsrf) = yalblw(j)  
          snow(i,nsrf) = ysnow(j)  
          qsurf(i,nsrf) = yqsurf(j)  
          rugos(i,nsrf) = yz0_new(j)  
          fluxlat(i,nsrf) = yfluxlat(j)  
 c$$$ pb         rugmer(i) = yrugm(j)  
          IF (nsrf .EQ. is_oce) then  
            rugmer(i) = yrugm(j)  
            rugos(i,nsrf) = yrugm(j)  
          endif    
 cIM cf JLD ??  
          agesno(i,nsrf) = yagesno(j)  
          fqcalving(i,nsrf) = y_fqcalving(j)          
          ffonte(i,nsrf) = y_ffonte(j)          
          cdragh(i) = cdragh(i) + ycoefh(j,1)  
          cdragm(i) = cdragm(i) + ycoefm(j,1)  
          dflux_t(i) = dflux_t(i) + y_dflux_t(j)  
          dflux_q(i) = dflux_q(i) + y_dflux_q(j)  
          zu1(i) = zu1(i) + yu1(j)  
          zv1(i) = zv1(i) + yv1(j)  
       END DO  
       IF ( nsrf .eq. is_ter ) THEN  
       DO j = 1, knon  
          i = ni(j)  
          qsol(i) = yqsol(j)  
       END DO  
       END IF  
       IF ( nsrf .eq. is_lic ) THEN  
         DO j = 1, knon  
           i = ni(j)  
           run_off_lic_0(i) = y_run_off_lic_0(j)  
         END DO  
       END IF  
 c$$$ PB ajout pour soil  
       ftsoil(:,:,nsrf) = 0.  
       DO k = 1, nsoilmx  
         DO j = 1, knon  
           i = ni(j)  
           ftsoil(i, k, nsrf) = ytsoil(j,k)  
         END DO  
       END DO  
 c  
       DO j = 1, knon  
       i = ni(j)  
       DO k = 1, klev  
          d_t(i,k) = d_t(i,k) + y_d_t(j,k)  
          d_q(i,k) = d_q(i,k) + y_d_q(j,k)  
 c$$$ PB        flux_t(i,k) = flux_t(i,k) + y_flux_t(j,k)  
 c$$$         flux_q(i,k) = flux_q(i,k) + y_flux_q(j,k)  
          d_u(i,k) = d_u(i,k) + y_d_u(j,k)  
          d_v(i,k) = d_v(i,k) + y_d_v(j,k)  
 c$$$  PB       flux_u(i,k) = flux_u(i,k) + y_flux_u(j,k)  
 c$$$         flux_v(i,k) = flux_v(i,k) + y_flux_v(j,k)  
          zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)  
       ENDDO  
       ENDDO  
 c  
 c  
 ccc diagnostic t,q a 2m et u, v a 10m  
 c  
       DO j=1, knon  
         i = ni(j)  
         uzon(j) = yu(j,1) + y_d_u(j,1)  
         vmer(j) = yv(j,1) + y_d_v(j,1)  
         tair1(j) = yt(j,1) + y_d_t(j,1)  
         qair1(j) = yq(j,1) + y_d_q(j,1)  
         zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1)))  
      &                   * (ypaprs(j,1)-ypplay(j,1))  
         tairsol(j) = yts(j) + y_d_ts(j)  
         rugo1(j) = yrugos(j)  
         IF(nsrf.EQ.is_oce) THEN  
          rugo1(j) = rugos(i,nsrf)  
         ENDIF  
         psfce(j)=ypaprs(j,1)  
         patm(j)=ypplay(j,1)  
 c  
         qairsol(j) = yqsurf(j)  
       ENDDO  
 c  
       CALL stdlevvar(klon, knon, nsrf, zxli,  
      &               uzon, vmer, tair1, qair1, zgeo1,  
      &               tairsol, qairsol, rugo1, psfce, patm,  
 cIM  &               yt2m, yq2m, yu10m)  
      &               yt2m, yq2m, yt10m, yq10m, yu10m, yustar)  
 cIM 081204 END  
 c  
 c  
       DO j=1, knon  
        i = ni(j)  
        t2m(i,nsrf)=yt2m(j)  
   
 c  
        q2m(i,nsrf)=yq2m(j)  
 c  
 c u10m, v10m : composantes du vent a 10m sans spirale de Ekman  
        u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)  
        v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)  
 c  
       ENDDO  
 c  
 cIM cf AM : pbl, HBTM  
       DO i = 1, knon  
          y_cd_h(i) = ycoefh(i,1)  
          y_cd_m(i) = ycoefm(i,1)  
       ENDDO  
 c     print*,'appel hbtm2'  
       CALL HBTM(knon, ypaprs, ypplay,  
      .          yt2m,yt10m,yq2m,yq10m,yustar,  
      .          y_flux_t,y_flux_q,yu,yv,yt,yq,  
      .          ypblh,ycapCL,yoliqCL,ycteiCL,ypblT,  
      .          ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)  
 c     print*,'fin hbtm2'  
 c  
       DO j=1, knon  
        i = ni(j)  
        pblh(i,nsrf)   = ypblh(j)  
        plcl(i,nsrf)   = ylcl(j)  
        capCL(i,nsrf)  = ycapCL(j)  
        oliqCL(i,nsrf) = yoliqCL(j)  
        cteiCL(i,nsrf) = ycteiCL(j)  
        pblT(i,nsrf)   = ypblT(j)  
        therm(i,nsrf)  = ytherm(j)  
        trmb1(i,nsrf)  = ytrmb1(j)  
        trmb2(i,nsrf)  = ytrmb2(j)  
        trmb3(i,nsrf)  = ytrmb3(j)  
       ENDDO  
 c  
   
       do j=1,knon  
          do k=1,klev+1  
          i=ni(j)  
          q2(i,k,nsrf)=yq2(j,k)  
          enddo  
       enddo  
 cIM "slab" ocean  
       IF(OCEAN.EQ.'slab  '.OR.OCEAN.EQ.'force ') THEN  
        IF (nsrf.EQ.is_oce) THEN  
         DO j = 1, knon  
 c on projette sur la grille globale  
          i = ni(j)  
          IF(pctsrf_new(i,is_oce).GT.epsfra) THEN  
           flux_o(i) = y_flux_o(j)  
          ELSE  
           flux_o(i) = 0.  
          ENDIF  
         ENDDO  
        ENDIF  
 c  
        IF (nsrf.EQ.is_sic) THEN  
         DO j = 1, knon  
          i = ni(j)  
 cIM 230604 on pondere lorsque l'on fait le bilan au sol :  flux_g(i) = y_flux_g(j)*ypct(j)  
          IF(pctsrf_new(i,is_sic).GT.epsfra) THEN  
           flux_g(i) = y_flux_g(j)  
          ELSE  
           flux_g(i) = 0.  
          ENDIF  
         ENDDO  
        ENDIF !nsrf.EQ.is_sic  
       ENDIF !OCEAN  
 c  
       IF(OCEAN.EQ.'slab  ') THEN  
        IF(nsrf.EQ.is_oce) then  
         tslab(1:klon) = ytslab(1:klon)  
         seaice(1:klon) = y_seaice(1:klon)  
        ENDIF !nsrf  
       ENDIF !OCEAN  
 99999 CONTINUE  
 C  
 C On utilise les nouvelles surfaces  
 C A rajouter: conservation de l'albedo  
 C  
       rugos(:,is_oce) = rugmer  
       pctsrf = pctsrf_new  
500    
501        RETURN            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
502        END         else
503              fsnow(:, nsrf) = 0.
504           end IF if_knon
505        END DO loop_surface
506    
507        ! On utilise les nouvelles surfaces
508        frugs(:, is_oce) = rugmer
509        pctsrf(:, is_oce) = pctsrf_new_oce
510        pctsrf(:, is_sic) = pctsrf_new_sic
511    
512        CALL histwrite_phy("run_off_lic", run_off_lic)
513        ftsol = ftsol + d_ts ! update surface temperature
514        CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
515        CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
516        CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
517        CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
518    
519      END SUBROUTINE pbl_surface
520    
521    end module pbl_surface_m

Legend:
Removed from v.3  
changed lines
  Added in v.332

  ViewVC Help
Powered by ViewVC 1.1.21