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

Diff of /trunk/phylmd/pbl_surface.f

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21