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

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

  ViewVC Help
Powered by ViewVC 1.1.21