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

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

  ViewVC Help
Powered by ViewVC 1.1.21