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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.21