/[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 13 by guez, Fri Jul 25 19:59:34 2008 UTC trunk/phylmd/Interface_surf/pbl_surface.f revision 308 by guez, Tue Sep 18 15:14:40 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, pplay, fsnow, qsurf, 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).      ! tableau des 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):: pplay(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):: qsurf(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(klon, nbsrf), 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, intent(in):: dtime      REAL therm(klon, nbsrf)
120        real date0      REAL plcl(klon, nbsrf)
121        integer, intent(in):: 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, intent(in):: 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)      REAL, intent(in):: tsol(:) ! (klon)
137        REAL seaice(klon), y_seaice(klon)  
138  cIM cf JLD      ! Local:
139        REAL y_fqcalving(klon), y_ffonte(klon)  
140        REAL fqcalving(klon,nbsrf), ffonte(klon,nbsrf)      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
141        REAL run_off_lic_0(klon), y_run_off_lic_0(klon)      REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
142    
143        REAL flux_u(klon,klev, nbsrf), flux_v(klon,klev, nbsrf)      ! la nouvelle repartition des surfaces sortie de l'interface
144        REAL rugmer(klon), agesno(klon,nbsrf)      REAL, save:: pctsrf_new_oce(klon)
145        real, intent(in):: rugoro(klon)      REAL, save:: pctsrf_new_sic(klon)
146        REAL cdragh(klon), cdragm(klon)  
147        integer jour            ! jour de l'annee en cours      REAL y_fqcalving(klon), y_ffonte(klon)
148        real rmu0(klon)         ! cosinus de l'angle solaire zenithal      real y_run_off_lic_0(klon), y_run_off_lic(klon)
149        REAL co2_ppm            ! taux CO2 atmosphere      REAL run_off_lic(klon) ! ruissellement total
150        LOGICAL, intent(in):: debut      REAL rugmer(klon)
151        logical, intent(in):: lafin      REAL ytsoil(klon, nsoilmx)
152        logical ok_veget      REAL yts(klon), ypct(klon), yz0_new(klon)
153        character(len=*), intent(IN):: ocean      real yrugos(klon) ! longueur de rugosite (en m)
154        integer npas, nexca      REAL yalb(klon)
155  c      REAL snow(klon), yqsurf(klon), yagesno(klon)
156        REAL pctsrf(klon,nbsrf)      real yqsol(klon) ! column-density of water in soil, in kg m-2
157        REAL ts(klon,nbsrf)      REAL yrain_fall(klon) ! liquid water mass flux (kg / m2 / s), positive down
158        REAL d_ts(klon,nbsrf)      REAL ysnow_fall(klon) ! solid water mass flux (kg / m2 / s), positive down
159        REAL snow(klon,nbsrf)      REAL yrugm(klon), radsol(klon), yrugoro(klon)
160        REAL qsurf(klon,nbsrf)      REAL yfluxlat(klon)
161        REAL evap(klon,nbsrf)      REAL y_d_ts(klon)
162        REAL albe(klon,nbsrf)      REAL y_d_t(klon, klev), y_d_q(klon, klev)
163        REAL alblw(klon,nbsrf)      REAL y_d_u(klon, klev), y_d_v(klon, klev)
164  c$$$ PB      REAL y_flux_t(klon), y_flux_q(klon)
165        REAL fluxlat(klon,nbsrf)      REAL y_flux_u(klon), y_flux_v(klon)
166  C      REAL y_dflux_t(klon), y_dflux_q(klon)
167        real rain_f(klon), snow_f(klon)      REAL ycoefh(klon, 2:klev), ycoefm(klon, 2:klev)
168        REAL fder(klon)      real ycdragh(klon), ycdragm(klon)
169  cIM cf. JLD   REAL sollw(klon), solsw(klon), sollwdown(klon)      REAL yu(klon, klev), yv(klon, klev)
170        REAL sollw(klon,nbsrf), solsw(klon,nbsrf), sollwdown(klon)      REAL yt(klon, klev), yq(klon, klev)
171        REAL rugos(klon,nbsrf)      REAL ypaprs(klon, klev + 1), ypplay(klon, klev), ydelp(klon, klev)
172  C la nouvelle repartition des surfaces sortie de l'interface      REAL yq2(klon, klev + 1)
173        REAL pctsrf_new(klon,nbsrf)      REAL delp(klon, klev)
174  cAA      INTEGER i, k, nsrf
175        REAL zcoefh(klon,klev)      INTEGER ni(klon), knon, j
176        REAL zu1(klon)  
177        REAL zv1(klon)      REAL pctsrf_pot(klon, nbsrf)
178  cAA      ! "pourcentage potentiel" pour tenir compte des \'eventuelles
179  c$$$ PB ajout pour soil      ! apparitions ou disparitions de la glace de mer
180        LOGICAL, intent(in):: soil_model  
181  cIM ajout seuils cdrm, cdrh      REAL yt2m(klon), yq2m(klon), wind10m(klon)
182        REAL cdmmax, cdhmax      REAL ustar(klon)
183  cIM: 261103  
184        REAL ksta, ksta_ter      REAL yt10m(klon), yq10m(klon)
185        LOGICAL ok_kzmin      REAL ypblh(klon)
186  cIM: 261103      REAL ylcl(klon)
187        REAL ftsoil(klon,nsoilmx,nbsrf)      REAL ycapcl(klon)
188        REAL ytsoil(klon,nsoilmx)      REAL yoliqcl(klon)
189        REAL qsol(klon)      REAL ycteicl(klon)
190  c======================================================================      REAL ypblt(klon)
191        EXTERNAL clqh, clvent, coefkz, calbeta, cltrac      REAL ytherm(klon)
192  c======================================================================      REAL u1(klon), v1(klon)
193        REAL yts(klon), yrugos(klon), ypct(klon), yz0_new(klon)      REAL tair1(klon), qair1(klon), tairsol(klon)
194        REAL yalb(klon)      REAL psfce(klon), patm(klon)
195        REAL yalblw(klon)      REAL zgeo1(klon)
196        REAL yu1(klon), yv1(klon)      REAL rugo1(klon)
197        real ysnow(klon), yqsurf(klon), yagesno(klon), yqsol(klon)      REAL zgeop(klon, klev)
198        real yrain_f(klon), ysnow_f(klon)  
199        real ysollw(klon), ysolsw(klon), ysollwdown(klon)      !------------------------------------------------------------
200        real yfder(klon), ytaux(klon), ytauy(klon)  
201        REAL yrugm(klon), yrads(klon),yrugoro(klon)      albsol = sum(falbe * pctsrf, dim = 2)
202  c$$$ PB  
203        REAL yfluxlat(klon)      ! R\'epartition sous maille des flux longwave et shortwave
204  C      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
205        REAL y_d_ts(klon)  
206        REAL y_d_t(klon, klev), y_d_q(klon, klev)      forall (nsrf = 1:nbsrf)
207        REAL y_d_u(klon, klev), y_d_v(klon, klev)         fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
208        REAL y_flux_t(klon,klev), y_flux_q(klon,klev)              * (tsol - ftsol(:, nsrf))
209        REAL y_flux_u(klon,klev), y_flux_v(klon,klev)         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
210        REAL y_dflux_t(klon), y_dflux_q(klon)      END forall
211        REAL ycoefh(klon,klev), ycoefm(klon,klev)  
212        REAL yu(klon,klev), yv(klon,klev)      ytherm = 0.
213        REAL yt(klon,klev), yq(klon,klev)  
214        REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)      DO k = 1, klev ! epaisseur de couche
215  c         DO i = 1, klon
216        LOGICAL ok_nonloc            delp(i, k) = paprs(i, k) - paprs(i, k + 1)
217        PARAMETER (ok_nonloc=.FALSE.)         END DO
218        REAL ycoefm0(klon,klev), ycoefh0(klon,klev)      END DO
219    
220  cIM 081204 hcl_Anne ? BEG      ! Initialization:
221        real yzlay(klon,klev),yzlev(klon,klev+1),yteta(klon,klev)      rugmer = 0.
222        real ykmm(klon,klev+1),ykmn(klon,klev+1)      cdragh = 0.
223        real ykmq(klon,klev+1)      cdragm = 0.
224        real yq2(klon,klev+1),q2(klon,klev+1,nbsrf)      dflux_t = 0.
225        real q2diag(klon,klev+1)      dflux_q = 0.
226  cIM 081204   real yustar(klon),y_cd_m(klon),y_cd_h(klon)      ypct = 0.
227  cIM 081204 hcl_Anne ? END      yrugos = 0.
228  c      ypaprs = 0.
229        REAL u1lay(klon), v1lay(klon)      ypplay = 0.
230        REAL delp(klon,klev)      ydelp = 0.
231        INTEGER i, k, nsrf      yrugoro = 0.
232  cAA   INTEGER it      d_ts = 0.
233        INTEGER ni(klon), knon, j      flux_t = 0.
234  c Introduction d'une variable "pourcentage potentiel" pour tenir compte      flux_q = 0.
235  c des eventuelles apparitions et/ou disparitions de la glace de mer      flux_u = 0.
236        REAL pctsrf_pot(klon,nbsrf)      flux_v = 0.
237              fluxlat = 0.
238  c======================================================================      d_t = 0.
239        REAL zx_alf1, zx_alf2 !valeur ambiante par extrapola.      d_q = 0.
240  c======================================================================      d_u = 0.
241  c      d_v = 0.
242  c maf pour sorties IOISPL en cas de debugagage      coefh = 0.
243  c      fqcalving = 0.
244        CHARACTER*80 cldebug      run_off_lic = 0.
245        SAVE cldebug  
246        CHARACTER*8 cl_surf(nbsrf)      ! Initialisation des "pourcentages potentiels". On consid\`ere ici qu'on
247        SAVE cl_surf      ! peut avoir potentiellement de la glace sur tout le domaine oc\'eanique
248        INTEGER nhoridbg, nidbg      ! (\`a affiner).
249        SAVE nhoridbg, nidbg  
250        INTEGER ndexbg(iim*(jjm+1))      pctsrf_pot(:, is_ter) = pctsrf(:, is_ter)
251        REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian      pctsrf_pot(:, is_lic) = pctsrf(:, is_lic)
252        REAL tabindx(klon)      pctsrf_pot(:, is_oce) = 1. - zmasq
253        REAL debugtab(iim,jjm+1)      pctsrf_pot(:, is_sic) = 1. - zmasq
254        LOGICAL first_appel  
255        SAVE first_appel      ! Tester si c'est le moment de lire le fichier:
256        DATA first_appel/.true./      if (mod(itap - 1, lmt_pas) == 0) then
257        LOGICAL:: debugindex = .false.         CALL interfoce_lim(julien, pctsrf_new_oce, pctsrf_new_sic)
258        integer idayref      endif
259        REAL t2m(klon,nbsrf), q2m(klon,nbsrf)  
260        REAL u10m(klon,nbsrf), v10m(klon,nbsrf)      ! Boucler sur toutes les sous-fractions du sol:
261  c  
262        REAL yt2m(klon), yq2m(klon), yu10m(klon)      loop_surface: DO nsrf = 1, nbsrf
263        REAL yustar(klon)         ! Define ni and knon:
264  c -- LOOP        
265         REAL yu10mx(klon)         ni = 0
266         REAL yu10my(klon)         knon = 0
267         REAL ywindsp(klon)  
268  c -- LOOP         DO i = 1, klon
269  c            ! Pour d\'eterminer le domaine \`a traiter, on utilise les surfaces
270        REAL yt10m(klon), yq10m(klon)            ! "potentielles"
271  cIM cf. AM : pbl, hbtm2 (Comme les autres diagnostics on cumule ds physic ce qui            IF (pctsrf_pot(i, nsrf) > epsfra) THEN
272  c   permet de sortir les grdeurs par sous surface)               knon = knon + 1
273        REAL pblh(klon,nbsrf)               ni(knon) = i
274        REAL plcl(klon,nbsrf)            END IF
275        REAL capCL(klon,nbsrf)         END DO
276        REAL oliqCL(klon,nbsrf)  
277        REAL cteiCL(klon,nbsrf)         if_knon: IF (knon /= 0) then
278        REAL pblT(klon,nbsrf)            DO j = 1, knon
279        REAL therm(klon,nbsrf)               i = ni(j)
280        REAL trmb1(klon,nbsrf)               ypct(j) = pctsrf(i, nsrf)
281        REAL trmb2(klon,nbsrf)               yts(j) = ftsol(i, nsrf)
282        REAL trmb3(klon,nbsrf)               snow(j) = fsnow(i, nsrf)
283        REAL ypblh(klon)               yqsurf(j) = qsurf(i, nsrf)
284        REAL ylcl(klon)               yalb(j) = falbe(i, nsrf)
285        REAL ycapCL(klon)               yrain_fall(j) = rain_fall(i)
286        REAL yoliqCL(klon)               ysnow_fall(j) = snow_fall(i)
287        REAL ycteiCL(klon)               yagesno(j) = agesno(i, nsrf)
288        REAL ypblT(klon)               yrugos(j) = frugs(i, nsrf)
289        REAL ytherm(klon)               yrugoro(j) = rugoro(i)
290        REAL ytrmb1(klon)               radsol(j) = fsolsw(i, nsrf) + fsollw(i, nsrf)
291        REAL ytrmb2(klon)               ypaprs(j, klev + 1) = paprs(i, klev + 1)
292        REAL ytrmb3(klon)               y_run_off_lic_0(j) = run_off_lic_0(i)
293        REAL y_cd_h(klon), y_cd_m(klon)            END DO
294  c     REAL ygamt(klon,2:klev) ! contre-gradient pour temperature  
295  c     REAL ygamq(klon,2:klev) ! contre-gradient pour humidite            ! For continent, copy soil water content
296        REAL uzon(klon), vmer(klon)            IF (nsrf == is_ter) yqsol(:knon) = qsol(ni(:knon))
297        REAL tair1(klon), qair1(klon), tairsol(klon)  
298        REAL psfce(klon), patm(klon)            ytsoil(:knon, :) = ftsoil(ni(:knon), :, nsrf)
299  c  
300        REAL qairsol(klon), zgeo1(klon)            DO k = 1, klev
301        REAL rugo1(klon)               DO j = 1, knon
302  c                  i = ni(j)
303        LOGICAL zxli ! utiliser un jeu de fonctions simples                  ypaprs(j, k) = paprs(i, k)
304        PARAMETER (zxli=.FALSE.)                  ypplay(j, k) = pplay(i, k)
305  c                  ydelp(j, k) = delp(i, k)
306        REAL zt, zqs, zdelta, zcor                  yu(j, k) = u(i, k)
307        REAL t_coup                  yv(j, k) = v(i, k)
308        PARAMETER(t_coup=273.15)                  yt(j, k) = t(i, k)
309  C                  yq(j, k) = q(i, k)
310        character (len = 20) :: modname = 'clmain'               END DO
       LOGICAL check  
       PARAMETER (check=.false.)  
   
   
 c initialisation Anne  
       ytherm(:) = 0.  
 C  
       if (check) THEN  
           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)  
           ENDDO  
           CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)  
           cldebug='sous_index'  
           CALL histbeg_totreg(cldebug, iim,zx_lon(:,1),jjm+1,  
      $        zx_lat(1,:),1,iim,1,jjm  
      $        +1, itau_phy,zjulian,dtime,nhoridbg,nidbg)  
 ! no vertical axis  
           cl_surf(1)='ter'  
           cl_surf(2)='lic'  
           cl_surf(3)='oce'  
           cl_surf(4)='sic'  
           DO nsrf=1,nbsrf  
             CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim,  
      $          jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)  
311            END DO            END DO
312            CALL histend(nidbg)  
313            CALL histsync(nidbg)            ! Calculer les géopotentiels de chaque couche:
314        ENDIF  
315              zgeop(:knon, 1) = RD * yt(:knon, 1) / (0.5 * (ypaprs(:knon, 1) &
316                   + ypplay(:knon, 1))) * (ypaprs(:knon, 1) - ypplay(:knon, 1))
317    
318              DO k = 2, klev
319                 zgeop(:knon, k) = zgeop(:knon, k - 1) + RD * 0.5 &
320                      * (yt(:knon, k - 1) + yt(:knon, k)) / ypaprs(:knon, k) &
321                      * (ypplay(:knon, k - 1) - ypplay(:knon, k))
322              ENDDO
323    
324              CALL cdrag(nsrf, sqrt(yu(:knon, 1)**2 + yv(:knon, 1)**2), &
325                   yt(:knon, 1), yq(:knon, 1), zgeop(:knon, 1), ypaprs(:knon, 1), &
326                   yts(:knon), yqsurf(:knon), yrugos(:knon), ycdragm(:knon), &
327                   ycdragh(:knon))
328    
329              IF (iflag_pbl == 1) THEN
330                 ycdragm(:knon) = max(ycdragm(:knon), 0.)
331                 ycdragh(:knon) = max(ycdragh(:knon), 0.)
332              end IF
333    
334              ! on met un seuil pour ycdragm et ycdragh
335              IF (nsrf == is_oce) THEN
336                 ycdragm(:knon) = min(ycdragm(:knon), cdmmax)
337                 ycdragh(:knon) = min(ycdragh(:knon), cdhmax)
338              END IF
339    
340              IF (iflag_pbl >= 6) yq2(:knon, :) = q2(ni(:knon), :, nsrf)
341              call coef_diff_turb(nsrf, ni(:knon), ypaprs(:knon, :), &
342                   ypplay(:knon, :), yu(:knon, :), yv(:knon, :), yq(:knon, :), &
343                   yt(:knon, :), yts(:knon), ycdragm(:knon), zgeop(:knon, :), &
344                   ycoefm(:knon, :), ycoefh(:knon, :), yq2(:knon, :))
345                        
346        DO k = 1, klev   ! epaisseur de couche            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
347        DO i = 1, klon                 ycdragm(:knon), yt(:knon, :), yu(:knon, :), ypaprs(:knon, :), &
348           delp(i,k) = paprs(i,k)-paprs(i,k+1)                 ypplay(:knon, :), ydelp(:knon, :), y_d_u(:knon, :), &
349        ENDDO                 y_flux_u(:knon))
350        ENDDO            CALL clvent(yu(:knon, 1), yv(:knon, 1), ycoefm(:knon, :), &
351        DO i = 1, klon  ! vent de la premiere couche                 ycdragm(:knon), yt(:knon, :), yv(:knon, :), ypaprs(:knon, :), &
352           zx_alf1 = 1.0                 ypplay(:knon, :), ydelp(:knon, :), y_d_v(:knon, :), &
353           zx_alf2 = 1.0 - zx_alf1                 y_flux_v(:knon))
354           u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2  
355           v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2            CALL clqh(julien, nsrf, ni(:knon), ytsoil(:knon, :), yqsol(:knon), &
356        ENDDO                 mu0(ni(:knon)), yrugos(:knon), yrugoro(:knon), yu(:knon, 1), &
357  c                 yv(:knon, 1), ycoefh(:knon, :), ycdragh(:knon), yt(:knon, :), &
358  c initialisation:                 yq(:knon, :), yts(:knon), ypaprs(:knon, :), ypplay(:knon, :), &
359  c                 ydelp(:knon, :), radsol(:knon), yalb(:knon), snow(:knon), &
360        DO i = 1, klon                 yqsurf(:knon), yrain_fall(:knon), ysnow_fall(:knon), &
361           rugmer(i) = 0.0                 yfluxlat(:knon), pctsrf_new_sic(ni(:knon)), yagesno(:knon), &
362           cdragh(i) = 0.0                 y_d_t(:knon, :), y_d_q(:knon, :), y_d_ts(:knon), &
363           cdragm(i) = 0.0                 yz0_new(:knon), y_flux_t(:knon), y_flux_q(:knon), &
364           dflux_t(i) = 0.0                 y_dflux_t(:knon), y_dflux_q(:knon), y_fqcalving(:knon), &
365           dflux_q(i) = 0.0                 y_ffonte(:knon), y_run_off_lic_0(:knon), y_run_off_lic(:knon))
366           zu1(i) = 0.0  
367           zv1(i) = 0.0            ! calculer la longueur de rugosite sur ocean
368        ENDDO  
369        ypct = 0.0            yrugm = 0.
370        yts = 0.0  
371        ysnow = 0.0            IF (nsrf == is_oce) THEN
372        yqsurf = 0.0               DO j = 1, knon
373        yalb = 0.0                  yrugm(j) = 0.018 * ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2) &
374        yalblw = 0.0                       / rg + 0.11 * 14E-6 &
375        yrain_f = 0.0                       / sqrt(ycdragm(j) * (yu(j, 1)**2 + yv(j, 1)**2))
376        ysnow_f = 0.0                  yrugm(j) = max(1.5E-05, yrugm(j))
377        yfder = 0.0               END DO
378        ytaux = 0.0            END IF
379        ytauy = 0.0  
380        ysolsw = 0.0            DO k = 1, klev
381        ysollw = 0.0               DO j = 1, knon
382        ysollwdown = 0.0                  i = ni(j)
383        yrugos = 0.0                  y_d_t(j, k) = y_d_t(j, k) * ypct(j)
384        yu1 = 0.0                  y_d_q(j, k) = y_d_q(j, k) * ypct(j)
385        yv1 = 0.0                  y_d_u(j, k) = y_d_u(j, k) * ypct(j)
386        yrads = 0.0                  y_d_v(j, k) = y_d_v(j, k) * ypct(j)
387        ypaprs = 0.0               END DO
388        ypplay = 0.0            END DO
389        ydelp = 0.0  
390        yu = 0.0            flux_t(ni(:knon), nsrf) = y_flux_t(:knon)
391        yv = 0.0            flux_q(ni(:knon), nsrf) = y_flux_q(:knon)
392        yt = 0.0            flux_u(ni(:knon), nsrf) = y_flux_u(:knon)
393        yq = 0.0            flux_v(ni(:knon), nsrf) = y_flux_v(:knon)
394        pctsrf_new = 0.0  
395        y_flux_u = 0.0            falbe(:, nsrf) = 0.
396        y_flux_v = 0.0            fsnow(:, nsrf) = 0.
397  C$$ PB            qsurf(:, nsrf) = 0.
398        y_dflux_t = 0.0            frugs(:, nsrf) = 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  
399            DO j = 1, knon            DO j = 1, knon
400              i = ni(j)               i = ni(j)
401              yqsol(j) = qsol(i)               d_ts(i, nsrf) = y_d_ts(j)
402                 falbe(i, nsrf) = yalb(j)
403                 fsnow(i, nsrf) = snow(j)
404                 qsurf(i, nsrf) = yqsurf(j)
405                 frugs(i, nsrf) = yz0_new(j)
406                 fluxlat(i, nsrf) = yfluxlat(j)
407                 IF (nsrf == is_oce) THEN
408                    rugmer(i) = yrugm(j)
409                    frugs(i, nsrf) = yrugm(j)
410                 END IF
411                 agesno(i, nsrf) = yagesno(j)
412                 fqcalving(i, nsrf) = y_fqcalving(j)
413                 ffonte(i, nsrf) = y_ffonte(j)
414                 cdragh(i) = cdragh(i) + ycdragh(j) * ypct(j)
415                 cdragm(i) = cdragm(i) + ycdragm(j) * ypct(j)
416                 dflux_t(i) = dflux_t(i) + y_dflux_t(j) * ypct(j)
417                 dflux_q(i) = dflux_q(i) + y_dflux_q(j) * ypct(j)
418            END DO            END DO
419        ELSE            IF (nsrf == is_ter) THEN
420            yqsol(:)=0.               qsol(ni(:knon)) = yqsol(:knon)
421        ENDIF            else IF (nsrf == is_lic) THEN
422  c$$$ PB ajour pour soil               DO j = 1, knon
423        DO k = 1, nsoilmx                  i = ni(j)
424          DO j = 1, knon                  run_off_lic_0(i) = y_run_off_lic_0(j)
425            i = ni(j)                  run_off_lic(i) = y_run_off_lic(j)
426            ytsoil(j,k) = ftsoil(i,k,nsrf)               END DO
427          END DO              END IF
428        END DO  
429        DO k = 1, klev            ftsoil(:, :, nsrf) = 0.
430        DO j = 1, knon            ftsoil(ni(:knon), :, nsrf) = ytsoil(:knon, :)
431        i = ni(j)  
432          ypaprs(j,k) = paprs(i,k)            DO j = 1, knon
433          ypplay(j,k) = pplay(i,k)               i = ni(j)
434          ydelp(j,k) = delp(i,k)               DO k = 1, klev
435          yu(j,k) = u(i,k)                  d_t(i, k) = d_t(i, k) + y_d_t(j, k)
436          yv(j,k) = v(i,k)                  d_q(i, k) = d_q(i, k) + y_d_q(j, k)
437          yt(j,k) = t(i,k)                  d_u(i, k) = d_u(i, k) + y_d_u(j, k)
438          yq(j,k) = q(i,k)                  d_v(i, k) = d_v(i, k) + y_d_v(j, k)
439        ENDDO               END DO
440        ENDDO            END DO
441  c  
442  c            forall (k = 2:klev) coefh(ni(:knon), k) &
443  c calculer Cdrag et les coefficients d'echange                 = coefh(ni(:knon), k) + ycoefh(:knon, k) * ypct(:knon)
444        CALL coefkz(nsrf, knon, ypaprs, ypplay,  
445  cIM 261103            ! diagnostic t, q a 2m et u, v a 10m
446       .     ksta, ksta_ter,  
447  cIM 261103            DO j = 1, knon
448       .            yts, yrugos, yu, yv, yt, yq,               i = ni(j)
449       .            yqsurf,               u1(j) = yu(j, 1) + y_d_u(j, 1)
450       .            ycoefm, ycoefh)               v1(j) = yv(j, 1) + y_d_v(j, 1)
451  cIM 081204 BEG               tair1(j) = yt(j, 1) + y_d_t(j, 1)
452  cCR test               qair1(j) = yq(j, 1) + y_d_q(j, 1)
453        if (iflag_pbl.eq.1) then               zgeo1(j) = rd * tair1(j) / (0.5 * (ypaprs(j, 1) + ypplay(j, &
454  cIM 081204 END                    1))) * (ypaprs(j, 1)-ypplay(j, 1))
455          CALL coefkz2(nsrf, knon, ypaprs, ypplay,yt,               tairsol(j) = yts(j) + y_d_ts(j)
456       .                  ycoefm0, ycoefh0)               rugo1(j) = yrugos(j)
457          DO k = 1, klev               IF (nsrf == is_oce) THEN
458          DO i = 1, knon                  rugo1(j) = frugs(i, nsrf)
459             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))               END IF
460             ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))               psfce(j) = ypaprs(j, 1)
461          ENDDO               patm(j) = ypplay(j, 1)
462          ENDDO            END DO
463        endif  
464  c            CALL stdlevvar(nsrf, u1(:knon), v1(:knon), tair1(:knon), qair1, &
465  cIM cf JLD : on seuille ycoefm et ycoefh                 zgeo1, tairsol, yqsurf(:knon), rugo1, psfce, patm, yt2m, yq2m, &
466        if (nsrf.eq.is_oce) then                 yt10m, yq10m, wind10m(:knon), ustar(:knon))
467           do j=1,knon  
468  c           ycoefm(j,1)=min(ycoefm(j,1),1.1E-3)            DO j = 1, knon
469              ycoefm(j,1)=min(ycoefm(j,1),cdmmax)               i = ni(j)
470  c           ycoefh(j,1)=min(ycoefh(j,1),1.1E-3)               t2m(i, nsrf) = yt2m(j)
471              ycoefh(j,1)=min(ycoefh(j,1),cdhmax)               q2m(i, nsrf) = yq2m(j)
472           enddo  
473        endif               u10m_srf(i, nsrf) = (wind10m(j) * u1(j)) &
474                      / sqrt(u1(j)**2 + v1(j)**2)
475  c               v10m_srf(i, nsrf) = (wind10m(j) * v1(j)) &
476  cIM: 261103                    / sqrt(u1(j)**2 + v1(j)**2)
477        if (ok_kzmin) THEN            END DO
478  cIM cf FH: 201103 BEG  
479  c   Calcul d'une diffusion minimale pour les conditions tres stables.            CALL hbtm(ypaprs, ypplay, yt2m, yq2m, ustar(:knon), y_flux_t(:knon), &
480        call coefkzmin(knon,ypaprs,ypplay,yu,yv,yt,yq,ycoefm                 y_flux_q(:knon), yu(:knon, :), yv(:knon, :), yt(:knon, :), &
481       .   ,ycoefm0,ycoefh0)                 yq(:knon, :), ypblh(:knon), ycapcl, yoliqcl, ycteicl, ypblt, &
482  c      call dump2d(iim,jjm-1,ycoefm(2:klon-1,2), 'KZ         ')                 ytherm, ylcl)
483  c      call dump2d(iim,jjm-1,ycoefm0(2:klon-1,2),'KZMIN      ')  
484              DO j = 1, knon
485         if ( 1.eq.1 ) then               i = ni(j)
486         DO k = 1, klev               pblh(i, nsrf) = ypblh(j)
487         DO i = 1, knon               plcl(i, nsrf) = ylcl(j)
488            ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))               capcl(i, nsrf) = ycapcl(j)
489            ycoefh(i,k) = MAX(ycoefh(i,k),ycoefh0(i,k))               oliqcl(i, nsrf) = yoliqcl(j)
490         ENDDO               cteicl(i, nsrf) = ycteicl(j)
491         ENDDO               pblt(i, nsrf) = ypblt(j)
492         endif               therm(i, nsrf) = ytherm(j)
493  cIM cf FH: 201103 END            END DO
494        endif !ok_kzmin  
495  cIM: 261103            IF (iflag_pbl >= 6) q2(ni(:knon), :, nsrf) = yq2(:knon, :)
496           else
497              fsnow(:, nsrf) = 0.
498        IF (iflag_pbl.ge.3) then         end IF if_knon
499        END DO loop_surface
500  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
501  c MELLOR ET YAMADA adapte a Mars Richard Fournier et Frederic Hourdin      ! On utilise les nouvelles surfaces
502  cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      frugs(:, is_oce) = rugmer
503        pctsrf(:, is_oce) = pctsrf_new_oce
504           yzlay(1:knon,1)=      pctsrf(:, is_sic) = pctsrf_new_sic
505       .   RD*yt(1:knon,1)/(0.5*(ypaprs(1:knon,1)+ypplay(1:knon,1)))  
506       .   *(ypaprs(1:knon,1)-ypplay(1:knon,1))/RG      CALL histwrite_phy("run_off_lic", run_off_lic)
507           do k=2,klev  
508              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  
           print *,'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 (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  
 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  
509    
510        RETURN  end module pbl_surface_m
       END  

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

  ViewVC Help
Powered by ViewVC 1.1.21