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

Legend:
Removed from v.12  
changed lines
  Added in v.305

  ViewVC Help
Powered by ViewVC 1.1.21