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

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

  ViewVC Help
Powered by ViewVC 1.1.21