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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21