source: LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90 @ 4674

Last change on this file since 4674 was 2886, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2865:2885 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 127.4 KB
Line 
1!
2! $Id: pbl_surface_mod.F90 2886 2017-05-20 07:41:16Z fhourdin $
3!
4MODULE pbl_surface_mod
5!
6! Planetary Boundary Layer and Surface module
7!
8! This module manages the calculation of turbulent diffusion in the boundary layer
9! and all interactions towards the differents sub-surfaces.
10!
11!
12  USE dimphy
13  USE mod_phys_lmdz_para,  ONLY : mpi_size
14  USE mod_grid_phy_lmdz,   ONLY : klon_glo
15  USE ioipsl
16  USE surface_data,        ONLY : type_ocean, ok_veget
17  USE surf_land_mod,       ONLY : surf_land
18  USE surf_landice_mod,    ONLY : surf_landice
19  USE surf_ocean_mod,      ONLY : surf_ocean
20  USE surf_seaice_mod,     ONLY : surf_seaice
21  USE cpl_mod,             ONLY : gath2cpl
22  USE climb_hq_mod,        ONLY : climb_hq_down, climb_hq_up
23  USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
24  USE coef_diff_turb_mod,  ONLY : coef_diff_turb
25
26
27  IMPLICIT NONE
28
29! Declaration of variables saved in restart file
30  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
31  !$OMP THREADPRIVATE(fder)
32  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: snow   ! snow at surface
33  !$OMP THREADPRIVATE(snow)
34  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface
35  !$OMP THREADPRIVATE(qsurf)
36  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature
37  !$OMP THREADPRIVATE(ftsoil)
38
39CONTAINS
40!
41!****************************************************************************************
42!
43  SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
44
45! This routine should be called after the restart file has been read.
46! This routine initialize the restart variables and does some validation tests
47! for the index of the different surfaces and tests the choice of type of ocean.
48
49    USE indice_sol_mod
50    USE print_control_mod, ONLY: lunout
51
52    INCLUDE "dimsoil.h"
53 
54! Input variables
55!****************************************************************************************
56    REAL, DIMENSION(klon), INTENT(IN)                 :: fder_rst
57    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: snow_rst
58    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: qsurf_rst
59    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
60
61 
62! Local variables
63!****************************************************************************************
64    INTEGER                       :: ierr
65    CHARACTER(len=80)             :: abort_message
66    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
67   
68
69!****************************************************************************************
70! Allocate and initialize module variables with fields read from restart file.
71!
72!****************************************************************************************   
73    ALLOCATE(fder(klon), stat=ierr)
74    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
75
76    ALLOCATE(snow(klon,nbsrf), stat=ierr)
77    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
78
79    ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
80    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
81
82    ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
83    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
84
85
86    fder(:)       = fder_rst(:)
87    snow(:,:)     = snow_rst(:,:)
88    qsurf(:,:)    = qsurf_rst(:,:)
89    ftsoil(:,:,:) = ftsoil_rst(:,:,:)
90
91
92!****************************************************************************************
93! Test for sub-surface indices
94!
95!****************************************************************************************
96    IF (is_ter /= 1) THEN
97      WRITE(lunout,*)" *** Warning ***"
98      WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
99      WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
100      abort_message="voir ci-dessus"
101      CALL abort_physic(modname,abort_message,1)
102    ENDIF
103
104    IF ( is_oce > is_sic ) THEN
105      WRITE(lunout,*)' *** Warning ***'
106      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
107      WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
108      WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
109      abort_message='voir ci-dessus'
110      CALL abort_physic(modname,abort_message,1)
111    ENDIF
112
113    IF ( is_lic > is_sic ) THEN
114      WRITE(lunout,*)' *** Warning ***'
115      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
116      WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
117      WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
118      abort_message='voir ci-dessus'
119      CALL abort_physic(modname,abort_message,1)
120    ENDIF
121
122!****************************************************************************************
123! Validation of ocean mode
124!
125!****************************************************************************************
126
127    IF (type_ocean /= 'slab  ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
128       WRITE(lunout,*)' *** Warning ***'
129       WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
130       abort_message='option pour l''ocean non valable'
131       CALL abort_physic(modname,abort_message,1)
132    ENDIF
133
134  END SUBROUTINE pbl_surface_init
135
136!****************************************************************************************
137
138
139  SUBROUTINE pbl_surface( &
140       dtime,     date0,     itap,     jour,          &
141       debut,     lafin,                              &
142       rlon,      rlat,      rugoro,   rmu0,          &
143       zsig,      lwdown_m,  pphi,     cldt,          &
144       rain_f,    snow_f,    solsw_m,  sollw_m,       &
145       gustiness,                                     &
146       t,         q,         u,        v,             &
147!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
148!!       t_x,       q_x,       t_w,      q_w,           &
149       wake_dlt,             wake_dlq,                &
150       wake_cstar,           wake_s,                  &
151!!!
152       pplay,     paprs,     pctsrf,                  &
153       ts,SFRWL,   alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
154       cdragh,    cdragm,   zu1,    zv1,              &
155       alb_dir_m,    alb_dif_m,  zxsens,   zxevap,    &
156       alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
157       zxtsol,    zxfluxlat, zt2m,     qsat2m,        &
158       d_t,       d_q,       d_u,      d_v, d_t_diss, &
159!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
160       d_t_w,     d_q_w,                              &
161       d_t_x,     d_q_x,                              &
162!!       d_wake_dlt,d_wake_dlq,                         &
163       zxsens_x,  zxfluxlat_x,zxsens_w,zxfluxlat_w,   &
164!!!
165!!! nrlmd le 13/06/2011
166       delta_tsurf,wake_dens,cdragh_x,cdragh_w,       &
167       cdragm_x,cdragm_w,kh,kh_x,kh_w,                &
168!!!
169       zcoefh,    zcoefm,    slab_wfbils,             &
170       qsol,    zq2m,      s_pblh,   s_plcl,        &
171!!!
172!!! jyg le 08/02/2012
173       s_pblh_x, s_plcl_x,   s_pblh_w, s_plcl_w,      &
174!!!
175       s_capCL,   s_oliqCL,  s_cteiCL, s_pblT,        &
176       s_therm,   s_trmb1,   s_trmb2,  s_trmb3,       &
177       zustar,zu10m,  zv10m,    fder_print,    &
178       zxqsurf,   rh2m,      zxfluxu,  zxfluxv,       &
179       z0m, z0h,   agesno,  sollw,    solsw,         &
180       d_ts,      evap,    fluxlat,  t2m,           &
181       wfbils,    wfbilo, wfevap, wfrain, wfsnow,   &
182       flux_t,   flux_u, flux_v,                    &
183       dflux_t,   dflux_q,   zxsnow,                  &
184!jyg<
185!!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
186       zxfluxt,   zxfluxq,   q2m,      flux_q, tke_x,   &
187!>jyg
188!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
189!!        tke_x,     tke_w                              &
190       wake_dltke                                     &
191!!!
192                        )
193!****************************************************************************************
194! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
195! Objet: interface de "couche limite" (diffusion verticale)
196!
197!AA REM:
198!AA-----
199!AA Tout ce qui a trait au traceurs est dans phytrac maintenant
200!AA pour l'instant le calcul de la couche limite pour les traceurs
201!AA se fait avec cltrac et ne tient pas compte de la differentiation
202!AA des sous-fraction de sol.
203!AA REM bis :
204!AA----------
205!AA Pour pouvoir extraire les coefficient d'echanges et le vent
206!AA dans la premiere couche, 3 champs supplementaires ont ete crees
207!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
208!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
209!AA si les informations des subsurfaces doivent etre prises en compte
210!AA il faudra sortir ces memes champs en leur ajoutant une dimension,
211!AA c'est a dire nbsrf (nbre de subsurface).
212!
213! Arguments:
214!
215! dtime----input-R- interval du temps (secondes)
216! itap-----input-I- numero du pas de temps
217! date0----input-R- jour initial
218! t--------input-R- temperature (K)
219! q--------input-R- vapeur d'eau (kg/kg)
220! u--------input-R- vitesse u
221! v--------input-R- vitesse v
222! wake_dlt-input-R- temperatre difference between (w) and (x) (K)
223! wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg)
224!wake_cstar-input-R- wake gust front speed (m/s)
225! wake_s---input-R- wake fractionnal area
226! ts-------input-R- temperature du sol (en Kelvin)
227! paprs----input-R- pression a intercouche (Pa)
228! pplay----input-R- pression au milieu de couche (Pa)
229! rlat-----input-R- latitude en degree
230! z0m, z0h ----input-R- longeur de rugosite (en m)
231! Martin
232! zsig-----input-R- slope
233! cldt-----input-R- total cloud fraction
234! pphi-----input-R- geopotentiel de chaque couche (g z) (reference sol)
235! Martin
236!
237! d_t------output-R- le changement pour "t"
238! d_q------output-R- le changement pour "q"
239! d_u------output-R- le changement pour "u"
240! d_v------output-R- le changement pour "v"
241! d_ts-----output-R- le changement pour "ts"
242! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
243!                    (orientation positive vers le bas)
244! tke_x---input/output-R- tke in the (x) region (kg/m**2/s)
245! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s)
246! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
247! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
248! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
249! dflux_t--output-R- derive du flux sensible
250! dflux_q--output-R- derive du flux latent
251! zu1------output-R- le vent dans la premiere couche
252! zv1------output-R- le vent dans la premiere couche
253! trmb1----output-R- deep_cape
254! trmb2----output-R- inhibition
255! trmb3----output-R- Point Omega
256! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
257! plcl-----output-R- Niveau de condensation
258! pblh-----output-R- HCL
259! pblT-----output-R- T au nveau HCL
260!
261    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
262    USE indice_sol_mod
263    USE time_phylmdz_mod, ONLY: day_ini,annee_ref,itau_phy
264    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
265    USE print_control_mod, ONLY: prt_level,lunout
266
267    IMPLICIT NONE
268
269    INCLUDE "dimsoil.h"
270    INCLUDE "YOMCST.h"
271    INCLUDE "YOETHF.h"
272    INCLUDE "FCTTRE.h"
273    INCLUDE "clesphys.h"
274    INCLUDE "compbl.h"
275    INCLUDE "flux_arp.h"
276!****************************************************************************************
277    REAL,                         INTENT(IN)        :: dtime   ! time interval (s)
278    REAL,                         INTENT(IN)        :: date0   ! initial day
279    INTEGER,                      INTENT(IN)        :: itap    ! time step
280    INTEGER,                      INTENT(IN)        :: jour    ! current day of the year
281    LOGICAL,                      INTENT(IN)        :: debut   ! true if first run step
282    LOGICAL,                      INTENT(IN)        :: lafin   ! true if last run step
283    REAL, DIMENSION(klon),        INTENT(IN)        :: rlon    ! longitudes in degrees
284    REAL, DIMENSION(klon),        INTENT(IN)        :: rlat    ! latitudes in degrees
285    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro  ! rugosity length
286    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosine of solar zenith angle
287    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
288    REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
289    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! net shortwave radiation at mean surface
290    REAL, DIMENSION(klon),        INTENT(IN)        :: sollw_m ! net longwave radiation at mean surface
291    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t       ! temperature (K)
292    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q       ! water vapour (kg/kg)
293    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u       ! u speed
294    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: v       ! v speed
295    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pplay   ! mid-layer pression (Pa)
296    REAL, DIMENSION(klon,klev+1), INTENT(IN)        :: paprs   ! pression between layers (Pa)
297    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf  ! sub-surface fraction
298! Martin
299    REAL, DIMENSION(klon),        INTENT(IN)        :: zsig    ! slope
300    REAL, DIMENSION(klon),        INTENT(IN)        :: lwdown_m ! downward longwave radiation at mean s   
301    REAL, DIMENSION(klon),        INTENT(IN)        :: gustiness ! gustiness
302
303    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud fraction
304    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pphi    ! geopotential (m2/s2)
305! Martin
306
307!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
308!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t_x       ! Temp\'erature hors poche froide
309!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t_w       ! Temp\'erature dans la poches froide
310!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q_x       !
311!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q_w       ! Pareil pour l'humidit\'e
312    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: wake_dlt  !temperature difference between (w) and (x) (K)
313    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: wake_dlq  !humidity difference between (w) and (x) (K)
314    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_s    ! Fraction de poches froides
315    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_cstar! Vitesse d'expansion des poches froides
316    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_dens
317!!!
318
319! Input/Output variables
320!****************************************************************************************
321    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
322    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: delta_tsurf !surface temperature difference between
323                                                                   !wake and off-wake regions
324!albedo SB >>>
325    REAL, DIMENSIOn(6),intent(in) :: SFRWL
326    REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT)     :: alb_dir,alb_dif
327!albedo SB <<<
328!jyg Pourquoi ustar et wstar sont-elles INOUT ?
329    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ustar   ! u* (m/s)
330    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)   :: wstar   ! w* (m/s)
331    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m    ! u speed at 10m
332    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: v10m    ! v speed at 10m
333!jyg<
334!!    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
335    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x
336!>jyg
337
338!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
339    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x
340!!!
341
342! Output variables
343!****************************************************************************************
344    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh     ! drag coefficient for T and Q
345    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm     ! drag coefficient for wind
346    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu1        ! u wind speed in first layer
347    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv1        ! v wind speed in first layer
348!albedo SB >>>
349    REAL, DIMENSION(klon, nsw),        INTENT(OUT)       :: alb_dir_m,alb_dif_m
350!albedo SB <<<
351    ! Martin
352    REAL, DIMENSION(klon),        INTENT(OUT)       :: alb3_lic
353    ! Martin
354    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens     ! sensible heat flux at surface with inversed sign
355                                                                  ! (=> positive sign upwards)
356    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxevap     ! water vapour flux at surface, positiv upwards
357    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol     ! temperature at surface, mean for each grid point
358!!! jyg le ???
359    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_t_w      !   !
360    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_q_w      !      !  Tendances dans les poches
361    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_t_x      !   !
362    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_q_x      !      !  Tendances hors des poches
363!!! jyg
364    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat  ! latent flux, mean for each grid point
365    REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m       ! temperature at 2m, mean for each grid point
366    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
367    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t        ! change in temperature
368    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t_diss       ! change in temperature
369    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_q        ! change in water vapour
370    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u        ! change in u speed
371    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_v        ! change in v speed
372
373    REAL, INTENT(OUT):: zcoefh(:, :, :) ! (klon, klev, nbsrf + 1)
374    ! coef for turbulent diffusion of T and Q, mean for each grid point
375
376    REAL, INTENT(OUT):: zcoefm(:, :, :) ! (klon, klev, nbsrf + 1)
377    ! coef for turbulent diffusion of U and V (?), mean for each grid point
378
379!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
380    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens_x   ! Flux sensible hors poche
381    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens_w   ! Flux sensible dans la poche
382    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat_x! Flux latent hors poche
383    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat_w! Flux latent dans la poche
384!!    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_wake_dlt
385!!    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_wake_dlq
386
387! Output only for diagnostics
388    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh_x
389    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh_w
390    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm_x
391    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm_w
392    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh
393    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh_x
394    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh_w
395!!!
396    REAL, DIMENSION(klon),        INTENT(OUT)       :: slab_wfbils! heat balance at surface only for slab at ocean points
397    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol     ! water height in the soil (mm)
398    REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m       ! water vapour at 2m, mean for each grid point
399    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh     ! height of the planetary boundary layer(HPBL)
400!!! jyg le 08/02/2012
401    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh_x   ! height of the PBL in the off-wake region
402    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh_w   ! height of the PBL in the wake region
403!!!
404    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl     ! condensation level
405!!! jyg le 08/02/2012
406    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl_x   ! condensation level in the off-wake region
407    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl_w   ! condensation level in the wake region
408!!!
409    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_capCL    ! CAPE of PBL
410    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_oliqCL   ! liquid water intergral of PBL
411    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_cteiCL   ! cloud top instab. crit. of PBL
412    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblT     ! temperature at PBLH
413    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_therm    ! thermal virtual temperature excess
414    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb1    ! deep cape, mean for each grid point
415    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb2    ! inhibition, mean for each grid point
416    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb3    ! point Omega, mean for each grid point
417    REAL, DIMENSION(klon),        INTENT(OUT)       :: zustar     ! u*
418    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu10m      ! u speed at 10m, mean for each grid point
419    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv10m      ! v speed at 10m, mean for each grid point
420    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
421    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf    ! humidity at surface, mean for each grid point
422    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m       ! relative humidity at 2m
423    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu    ! u wind tension, mean for each grid point
424    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxv    ! v wind tension, mean for each grid point
425    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)     :: z0m,z0h      ! rugosity length (m)
426    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: agesno   ! age of snow at surface
427    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: solsw      ! net shortwave radiation at surface
428    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw      ! net longwave radiation at surface
429    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts       ! change in temperature at surface
430    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)       :: evap     ! evaporation at surface
431    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat    ! latent flux
432    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m        ! temperature at 2 meter height
433    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbils     ! heat balance at surface
434    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbilo     ! water balance at surface weighted by srf
435    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfevap     ! water balance (evap) at surface weighted by srf
436    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfrain     ! water balance (rain) at surface weighted by srf
437    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfsnow     ! water balance (snow) at surface weighted by srf
438    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t     ! sensible heat flux (CpT) J/m**2/s (W/m**2)
439                                                                  ! positve orientation downwards
440    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u     ! u wind tension (kg m/s)/(m**2 s) or Pascal
441    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v     ! v wind tension (kg m/s)/(m**2 s) or Pascal
442
443! Output not needed
444    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_t    ! change of sensible heat flux
445    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_q    ! change of water vapour flux
446    REAL, DIMENSION(klon),       INTENT(OUT)        :: zxsnow     ! snow at surface, mean for each grid point
447    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt    ! sensible heat flux, mean for each grid point
448    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxq    ! water vapour flux, mean for each grid point
449    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m        ! water vapour at 2 meter height
450    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
451
452! Martin
453! sisvat
454    REAL, DIMENSION(klon),       INTENT(OUT)        :: qsnow      ! snow water content
455    REAL, DIMENSION(klon),       INTENT(OUT)        :: snowhgt    ! snow height
456    REAL, DIMENSION(klon),       INTENT(OUT)        :: to_ice     ! snow passed to ice
457    REAL, DIMENSION(klon),       INTENT(OUT)        :: sissnow    ! snow in snow model
458    REAL, DIMENSION(klon),       INTENT(OUT)        :: runoff     ! runoff on land ice
459! Martin
460
461! Local variables with attribute SAVE
462!****************************************************************************************
463    INTEGER, SAVE                            :: nhoridbg, nidbg   ! variables for IOIPSL
464!$OMP THREADPRIVATE(nhoridbg, nidbg)
465    LOGICAL, SAVE                            :: debugindex=.FALSE.
466!$OMP THREADPRIVATE(debugindex)
467    LOGICAL, SAVE                            :: first_call=.TRUE.
468!$OMP THREADPRIVATE(first_call)
469    CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
470!$OMP THREADPRIVATE(cl_surf)
471
472! Other local variables
473!****************************************************************************************
474    INTEGER                            :: iflag_split
475    INTEGER                            :: i, k, nsrf
476    INTEGER                            :: knon, j
477    INTEGER                            :: idayref
478    INTEGER , DIMENSION(klon)          :: ni
479    REAL                               :: yt1_new
480    REAL                               :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
481    REAL                               :: amn, amx
482    REAL                               :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
483    REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
484    REAL, DIMENSION(klon)              :: yts, yz0m, yz0h, ypct
485!albedo SB >>>
486    REAL, DIMENSION(klon)              :: yalb,yalb_vis
487!albedo SB <<<
488    REAL, DIMENSION(klon)              :: yu1, yv1
489    REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
490    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f
491    REAL, DIMENSION(klon)              :: ysolsw, ysollw
492    REAL, DIMENSION(klon)              :: yfder
493    REAL, DIMENSION(klon)              :: yrugoro
494    REAL, DIMENSION(klon)              :: yfluxlat
495    REAL, DIMENSION(klon)              :: y_d_ts
496    REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
497    REAL, DIMENSION(klon)              :: y_dflux_t, y_dflux_q
498    REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
499    REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
500    REAL, DIMENSION(klon)              :: yustar
501    REAL, DIMENSION(klon)              :: ywstar
502    REAL, DIMENSION(klon)              :: ywindsp
503    REAL, DIMENSION(klon)              :: yt10m, yq10m
504    REAL, DIMENSION(klon)              :: ypblh
505    REAL, DIMENSION(klon)              :: ylcl
506    REAL, DIMENSION(klon)              :: ycapCL
507    REAL, DIMENSION(klon)              :: yoliqCL
508    REAL, DIMENSION(klon)              :: ycteiCL
509    REAL, DIMENSION(klon)              :: ypblT
510    REAL, DIMENSION(klon)              :: ytherm
511    REAL, DIMENSION(klon)              :: ytrmb1
512    REAL, DIMENSION(klon)              :: ytrmb2
513    REAL, DIMENSION(klon)              :: ytrmb3
514    REAL, DIMENSION(klon)              :: uzon, vmer
515    REAL, DIMENSION(klon)              :: tair1, qair1, tairsol
516    REAL, DIMENSION(klon)              :: psfce, patm
517    REAL, DIMENSION(klon)              :: qairsol, zgeo1, speed, zri1, pref !speed, zri1, pref, added by Fuxing WANG, 04/03/2015
518    REAL, DIMENSION(klon)              :: rugo1
519    REAL, DIMENSION(klon)              :: yfluxsens
520    REAL, DIMENSION(klon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
521    REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
522    REAL, DIMENSION(klon)              :: ypsref
523    REAL, DIMENSION(klon)              :: yevap, ytsurf_new, yalb3_new
524!albedo SB >>>
525    REAL, DIMENSION(klon,nsw)          :: yalb_dir_new, yalb_dif_new
526!albedo SB <<<
527    REAL, DIMENSION(klon)              :: ztsol
528    REAL, DIMENSION(klon)              :: meansqT ! mean square deviation of subsurface temperatures
529    REAL, DIMENSION(klon)              :: alb_m  ! mean albedo for whole SW interval
530    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q, y_d_t_diss
531    REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
532    REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q
533    REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
534    REAL, DIMENSION(klon,klev)         :: ycoefh, ycoefm,ycoefq
535    REAL, DIMENSION(klon)              :: ycdragh, ycdragm
536    REAL, DIMENSION(klon,klev)         :: yu, yv
537    REAL, DIMENSION(klon,klev)         :: yt, yq
538    REAL, DIMENSION(klon,klev)         :: ypplay, ydelp
539    REAL, DIMENSION(klon,klev)         :: delp
540    REAL, DIMENSION(klon,klev+1)       :: ypaprs
541    REAL, DIMENSION(klon,klev+1)       :: ytke
542    REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
543    CHARACTER(len=80)                  :: abort_message
544    CHARACTER(len=20)                  :: modname = 'pbl_surface'
545    LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
546    LOGICAL, PARAMETER                 :: check=.FALSE.
547
548!!! nrlmd le 02/05/2011
549!!! jyg le 07/02/2012
550    REAL, DIMENSION(klon)              :: ywake_s, ywake_cstar, ywake_dens
551!!!
552    REAL, DIMENSION(klon,klev+1)       :: ytke_x, ytke_w
553    REAL, DIMENSION(klon,klev+1)       :: ywake_dltke
554    REAL, DIMENSION(klon,klev)         :: yu_x, yv_x, yu_w, yv_w
555    REAL, DIMENSION(klon,klev)         :: yt_x, yq_x, yt_w, yq_w
556    REAL, DIMENSION(klon,klev)         :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w
557    REAL, DIMENSION(klon,klev)         :: ycoefq_x, ycoefq_w
558    REAL, DIMENSION(klon)              :: ycdragh_x, ycdragm_x, ycdragh_w, ycdragm_w
559    REAL, DIMENSION(klon)              :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x
560    REAL, DIMENSION(klon)              :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w
561    REAL, DIMENSION(klon)              :: AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x
562    REAL, DIMENSION(klon)              :: AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w
563    REAL, DIMENSION(klon)              :: y_flux_t1_x, y_flux_q1_x, y_flux_t1_w, y_flux_q1_w
564    REAL, DIMENSION(klon)              :: y_flux_u1_x, y_flux_v1_x, y_flux_u1_w, y_flux_v1_w
565    REAL, DIMENSION(klon,klev)         :: y_flux_t_x, y_flux_q_x, y_flux_t_w, y_flux_q_w
566    REAL, DIMENSION(klon,klev)         :: y_flux_u_x, y_flux_v_x, y_flux_u_w, y_flux_v_w
567    REAL, DIMENSION(klon)              :: yfluxlat_x, yfluxlat_w
568    REAL, DIMENSION(klon,klev)         :: y_d_t_x, y_d_q_x, y_d_t_w, y_d_q_w
569    REAL, DIMENSION(klon,klev)         :: y_d_t_diss_x, y_d_t_diss_w
570    REAL, DIMENSION(klon,klev)         :: d_t_diss_x, d_t_diss_w
571    REAL, DIMENSION(klon,klev)         :: y_d_u_x, y_d_v_x, y_d_u_w, y_d_v_w
572    REAL, DIMENSION(klon, klev, nbsrf) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w
573    REAL, DIMENSION(klon, klev, nbsrf) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w
574    REAL, DIMENSION(klon, nbsrf)       :: fluxlat_x, fluxlat_w
575    REAL, DIMENSION(klon, klev)        :: zxfluxt_x, zxfluxq_x, zxfluxt_w, zxfluxq_w
576    REAL, DIMENSION(klon, klev)        :: zxfluxu_x, zxfluxv_x, zxfluxu_w, zxfluxv_w
577    REAL                               :: zx_qs_surf, zcor_surf, zdelta_surf
578    REAL, DIMENSION(klon)              :: ytsurf_th, yqsatsurf
579    REAL, DIMENSION(klon)              :: ybeta
580    REAL, DIMENSION(klon, klev)        :: d_u_x
581    REAL, DIMENSION(klon, klev)        :: d_u_w
582    REAL, DIMENSION(klon, klev)        :: d_v_x
583    REAL, DIMENSION(klon, klev)        :: d_v_w
584
585    REAL, DIMENSION(klon,klev)         :: CcoefH, CcoefQ, DcoefH, DcoefQ
586    REAL, DIMENSION(klon,klev)         :: CcoefU, CcoefV, DcoefU, DcoefV
587    REAL, DIMENSION(klon,klev)         :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x
588    REAL, DIMENSION(klon,klev)         :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w
589    REAL, DIMENSION(klon,klev)         :: CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x
590    REAL, DIMENSION(klon,klev)         :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w
591    REAL, DIMENSION(klon,klev)         :: Kcoef_hq, Kcoef_m, gama_h, gama_q
592    REAL, DIMENSION(klon,klev)         :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x
593    REAL, DIMENSION(klon,klev)         :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
594    REAL, DIMENSION(klon)              :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w
595!!!
596!!!jyg le 08/02/2012
597    REAL, DIMENSION(klon, nbsrf)       :: windsp
598!
599    REAL, DIMENSION(klon, nbsrf)       :: t2m_x
600    REAL, DIMENSION(klon, nbsrf)       :: q2m_x
601    REAL, DIMENSION(klon)              :: rh2m_x
602    REAL, DIMENSION(klon)              :: qsat2m_x
603    REAL, DIMENSION(klon, nbsrf)       :: u10m_x
604    REAL, DIMENSION(klon, nbsrf)       :: v10m_x
605    REAL, DIMENSION(klon, nbsrf)       :: ustar_x
606    REAL, DIMENSION(klon, nbsrf)       :: wstar_x
607!             
608    REAL, DIMENSION(klon, nbsrf)       :: pblh_x
609    REAL, DIMENSION(klon, nbsrf)       :: plcl_x
610    REAL, DIMENSION(klon, nbsrf)       :: capCL_x
611    REAL, DIMENSION(klon, nbsrf)       :: oliqCL_x
612    REAL, DIMENSION(klon, nbsrf)       :: cteiCL_x
613    REAL, DIMENSION(klon, nbsrf)       :: pblt_x
614    REAL, DIMENSION(klon, nbsrf)       :: therm_x
615    REAL, DIMENSION(klon, nbsrf)       :: trmb1_x
616    REAL, DIMENSION(klon, nbsrf)       :: trmb2_x
617    REAL, DIMENSION(klon, nbsrf)       :: trmb3_x
618!
619    REAL, DIMENSION(klon, nbsrf)       :: t2m_w
620    REAL, DIMENSION(klon, nbsrf)       :: q2m_w
621    REAL, DIMENSION(klon)              :: rh2m_w
622    REAL, DIMENSION(klon)              :: qsat2m_w
623    REAL, DIMENSION(klon, nbsrf)       :: u10m_w
624    REAL, DIMENSION(klon, nbsrf)       :: v10m_w
625    REAL, DIMENSION(klon, nbsrf)       :: ustar_w
626    REAL, DIMENSION(klon, nbsrf)       :: wstar_w
627!                           
628    REAL, DIMENSION(klon, nbsrf)       :: pblh_w
629    REAL, DIMENSION(klon, nbsrf)       :: plcl_w
630    REAL, DIMENSION(klon, nbsrf)       :: capCL_w
631    REAL, DIMENSION(klon, nbsrf)       :: oliqCL_w
632    REAL, DIMENSION(klon, nbsrf)       :: cteiCL_w
633    REAL, DIMENSION(klon, nbsrf)       :: pblt_w
634    REAL, DIMENSION(klon, nbsrf)       :: therm_w
635    REAL, DIMENSION(klon, nbsrf)       :: trmb1_w
636    REAL, DIMENSION(klon, nbsrf)       :: trmb2_w
637    REAL, DIMENSION(klon, nbsrf)       :: trmb3_w
638!
639    REAL, DIMENSION(klon)       :: yt2m_x
640    REAL, DIMENSION(klon)       :: yq2m_x
641    REAL, DIMENSION(klon)       :: yt10m_x
642    REAL, DIMENSION(klon)       :: yq10m_x
643    REAL, DIMENSION(klon)       :: yu10m_x
644    REAL, DIMENSION(klon)       :: yv10m_x
645    REAL, DIMENSION(klon)       :: yustar_x
646    REAL, DIMENSION(klon)       :: ywstar_x
647!             
648    REAL, DIMENSION(klon)       :: ypblh_x
649    REAL, DIMENSION(klon)       :: ylcl_x
650    REAL, DIMENSION(klon)       :: ycapCL_x
651    REAL, DIMENSION(klon)       :: yoliqCL_x
652    REAL, DIMENSION(klon)       :: ycteiCL_x
653    REAL, DIMENSION(klon)       :: ypblt_x
654    REAL, DIMENSION(klon)       :: ytherm_x
655    REAL, DIMENSION(klon)       :: ytrmb1_x
656    REAL, DIMENSION(klon)       :: ytrmb2_x
657    REAL, DIMENSION(klon)       :: ytrmb3_x
658!
659    REAL, DIMENSION(klon)       :: yt2m_w
660    REAL, DIMENSION(klon)       :: yq2m_w
661    REAL, DIMENSION(klon)       :: yt10m_w
662    REAL, DIMENSION(klon)       :: yq10m_w
663    REAL, DIMENSION(klon)       :: yu10m_w
664    REAL, DIMENSION(klon)       :: yv10m_w
665    REAL, DIMENSION(klon)       :: yustar_w
666    REAL, DIMENSION(klon)       :: ywstar_w
667!                       
668    REAL, DIMENSION(klon)       :: ypblh_w
669    REAL, DIMENSION(klon)       :: ylcl_w
670    REAL, DIMENSION(klon)       :: ycapCL_w
671    REAL, DIMENSION(klon)       :: yoliqCL_w
672    REAL, DIMENSION(klon)       :: ycteiCL_w
673    REAL, DIMENSION(klon)       :: ypblt_w
674    REAL, DIMENSION(klon)       :: ytherm_w
675    REAL, DIMENSION(klon)       :: ytrmb1_w
676    REAL, DIMENSION(klon)       :: ytrmb2_w
677    REAL, DIMENSION(klon)       :: ytrmb3_w
678!
679    REAL, DIMENSION(klon)              :: uzon_x, vmer_x, speed_x, zri1_x, pref_x !speed_x, zri1_x, pref_x, added by Fuxing WANG, 04/03/2015
680    REAL, DIMENSION(klon)              :: zgeo1_x, tair1_x, qair1_x, tairsol_x
681!
682    REAL, DIMENSION(klon)              :: uzon_w, vmer_w, speed_w, zri1_w, pref_w !speed_w, zri1_w, pref_w, added by Fuxing WANG, 04/03/2015
683    REAL, DIMENSION(klon)              :: zgeo1_w, tair1_w, qair1_w, tairsol_w
684
685!!! jyg le 25/03/2013
686!!    Variables intermediaires pour le raccord des deux colonnes \`a la surface
687    REAL   ::   dd_Ch
688    REAL   ::   dd_Cm
689    REAL   ::   dd_Kh
690    REAL   ::   dd_Km
691    REAL   ::   dd_u
692    REAL   ::   dd_v
693    REAL   ::   dd_t
694    REAL   ::   dd_q
695    REAL   ::   dd_AH
696    REAL   ::   dd_AQ
697    REAL   ::   dd_AU
698    REAL   ::   dd_AV
699    REAL   ::   dd_BH
700    REAL   ::   dd_BQ
701    REAL   ::   dd_BU
702    REAL   ::   dd_BV
703
704    REAL   ::   dd_KHp
705    REAL   ::   dd_KQp
706    REAL   ::   dd_KUp
707    REAL   ::   dd_KVp
708
709!!!
710!!! nrlmd le 13/06/2011
711    REAL, DIMENSION(klon)              :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1
712    REAL, DIMENSION(klon)              :: y_delta_tsurf,delta_coef,tau_eq
713    REAL, PARAMETER                    :: facteur=2./sqrt(3.14)
714    REAL, PARAMETER                    :: effusivity=2000.
715    REAL, DIMENSION(klon)              :: ytsurf_th_x,ytsurf_th_w,yqsatsurf_x,yqsatsurf_w
716    REAL, DIMENSION(klon)              :: ydtsurf_th
717    REAL                               :: zdelta_surf_x,zdelta_surf_w,zx_qs_surf_x,zx_qs_surf_w
718    REAL                               :: zcor_surf_x,zcor_surf_w
719    REAL                               :: mod_wind_x, mod_wind_w
720    REAL                               :: rho1
721    REAL, DIMENSION(klon)              :: Kech_h           ! Coefficient d'echange pour l'energie
722    REAL, DIMENSION(klon)              :: Kech_h_x, Kech_h_w
723    REAL, DIMENSION(klon)              :: Kech_m
724    REAL, DIMENSION(klon)              :: Kech_m_x, Kech_m_w
725    REAL, DIMENSION(klon)              :: yts_x,yts_w
726    REAL, DIMENSION(klon)              :: Kech_Hp, Kech_H_xp, Kech_H_wp
727    REAL, DIMENSION(klon)              :: Kech_Qp, Kech_Q_xp, Kech_Q_wp
728    REAL, DIMENSION(klon)              :: Kech_Up, Kech_U_xp, Kech_U_wp
729    REAL, DIMENSION(klon)              :: Kech_Vp, Kech_V_xp, Kech_V_wp
730
731    REAL                               :: vent
732
733
734
735
736!!!
737
738! For debugging with IOIPSL
739    INTEGER, DIMENSION(nbp_lon*nbp_lat)    :: ndexbg
740    REAL                               :: zjulian
741    REAL, DIMENSION(klon)              :: tabindx
742    REAL, DIMENSION(nbp_lon,nbp_lat)         :: zx_lon, zx_lat
743    REAL, DIMENSION(nbp_lon,nbp_lat)         :: debugtab
744
745
746    REAL, DIMENSION(klon,nbsrf)        :: pblh         ! height of the planetary boundary layer
747    REAL, DIMENSION(klon,nbsrf)        :: plcl         ! condensation level
748    REAL, DIMENSION(klon,nbsrf)        :: capCL
749    REAL, DIMENSION(klon,nbsrf)        :: oliqCL
750    REAL, DIMENSION(klon,nbsrf)        :: cteiCL
751    REAL, DIMENSION(klon,nbsrf)        :: pblT
752    REAL, DIMENSION(klon,nbsrf)        :: therm
753    REAL, DIMENSION(klon,nbsrf)        :: trmb1        ! deep cape
754    REAL, DIMENSION(klon,nbsrf)        :: trmb2        ! inhibition
755    REAL, DIMENSION(klon,nbsrf)        :: trmb3        ! point Omega
756    REAL, DIMENSION(klon,nbsrf)        :: zx_rh2m, zx_qsat2m
757    REAL, DIMENSION(klon,nbsrf)        :: zx_t1
758    REAL, DIMENSION(klon, nbsrf)       :: alb          ! mean albedo for whole SW interval
759    REAL, DIMENSION(klon)              :: ylwdown      ! jg : temporary (ysollwdown)
760    REAL, DIMENSION(klon)              :: ygustiness      ! jg : temporary (ysollwdown)
761
762    REAL                               :: zx_qs1, zcor1, zdelta1
763
764    ! Martin
765    REAL, DIMENSION(klon, nbsrf)       :: sollwd ! net longwave radiation at surface
766    REAL, DIMENSION(klon)              :: ytoice
767    REAL, DIMENSION(klon)              :: ysnowhgt, yqsnow, ysissnow, yrunoff
768    REAL, DIMENSION(klon)              :: yzsig
769    REAL, DIMENSION(klon,klev)         :: ypphi
770    REAL, DIMENSION(klon)              :: ycldt
771    REAL, DIMENSION(klon)              :: yrmu0
772    ! Martin
773
774!****************************************************************************************
775! End of declarations
776!****************************************************************************************
777
778      IF (prt_level >=10) print *,' -> pbl_surface, itap ',itap
779!
780!!jyg      iflag_split = mod(iflag_pbl_split,2)
781      iflag_split = mod(iflag_pbl_split,10)
782
783!****************************************************************************************
784! 1) Initialisation and validation tests
785!    Only done first time entering this subroutine
786!
787!****************************************************************************************
788
789    IF (first_call) THEN
790       print*,'PBL SURFACE AVEC GUSTINESS'
791       first_call=.FALSE.
792     
793       ! Initialize ok_flux_surf (for 1D model)
794       if (klon_glo>1) ok_flux_surf=.FALSE.
795       
796       ! Initilize debug IO
797       IF (debugindex .AND. mpi_size==1) THEN
798          ! initialize IOIPSL output
799          idayref = day_ini
800          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
801          CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlon,zx_lon)
802          DO i = 1, nbp_lon
803             zx_lon(i,1) = rlon(i+1)
804             zx_lon(i,nbp_lat) = rlon(i+1)
805          ENDDO
806          CALL gr_fi_ecrit(1,klon,nbp_lon,nbp_lat,rlat,zx_lat)
807          CALL histbeg("sous_index",nbp_lon,zx_lon(:,1),nbp_lat,zx_lat(1,:), &
808               1,nbp_lon,1,nbp_lat, &
809               itau_phy,zjulian,dtime,nhoridbg,nidbg)
810          ! no vertical axis
811          cl_surf(1)='ter'
812          cl_surf(2)='lic'
813          cl_surf(3)='oce'
814          cl_surf(4)='sic'
815          DO nsrf=1,nbsrf
816             CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",nbp_lon, &
817                  nbp_lat,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
818          END DO
819
820          CALL histend(nidbg)
821          CALL histsync(nidbg)
822
823       END IF
824       
825    ENDIF
826         
827!****************************************************************************************
828! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
829! instead of ORCHIDEE)
830    IF (qsol0>=0.) THEN
831      PRINT*,'WARNING : On impose qsol=',qsol0
832      qsol(:)=qsol0
833    ENDIF
834!****************************************************************************************
835
836!****************************************************************************************
837! 2) Initialization to zero
838!****************************************************************************************
839!
840! 2a) Initialization of all argument variables with INTENT(OUT)
841!****************************************************************************************
842 cdragh(:)=0. ; cdragm(:)=0.
843 zu1(:)=0. ; zv1(:)=0.
844!albedo SB >>>
845  alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.
846!albedo SB <<<
847 zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0.
848 d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
849 zxfluxlat(:)=0.
850 zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
851 d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
852 zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
853 cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
854 kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.
855 slab_wfbils(:)=0.
856 s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.
857 s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.
858 s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.
859 s_therm(:)=0.
860 s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.
861 zustar(:)=0.
862 zu10m(:)=0. ; zv10m(:)=0.
863 fder_print(:)=0.
864 zxqsurf(:)=0.
865 zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
866 solsw(:,:)=0. ; sollw(:,:)=0.
867 d_ts(:,:)=0.
868 evap(:,:)=0.
869 fluxlat(:,:)=0.
870 wfbils(:,:)=0. ; wfbilo(:,:)=0.
871 wfevap(:,:)=0. ; wfrain(:,:)=0. ; wfsnow(:,:)=0.
872 flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
873 dflux_t(:)=0. ; dflux_q(:)=0.
874 zxsnow(:)=0.
875 zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.
876 qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
877 runoff(:)=0.
878    IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
879       zcoefh(:,:,:) = 0.0
880       zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
881       zcoefm(:,:,:) = 0.0
882       zcoefm(:,1,:) = 999999. !
883    ELSE
884      zcoefm(:,:,is_ave)=0.
885      zcoefh(:,:,is_ave)=0.
886    ENDIF
887!!
888!  The components "is_ave" of tke_x and wake_deltke are "OUT" variables
889!jyg<
890!!    tke(:,:,is_ave)=0.
891    tke_x(:,:,is_ave)=0.
892    wake_dltke(:,:,is_ave)=0.
893!>jyg
894!!! jyg le 23/02/2013
895    t2m(:,:)       = 999999.     ! t2m and q2m are meaningfull only over sub-surfaces
896    q2m(:,:)       = 999999.     ! actually present in the grid cell.
897!!!
898    rh2m(:) = 0. ; qsat2m(:) = 0.
899!!!
900!!! jyg le 10/02/2012
901    rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.
902!!!
903
904! 2b) Initialization of all local variables that will be compressed later
905!****************************************************************************************
906!!    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
907    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0
908!!    zv1 = 0.0     ; yqsurf = 0.0
909!albedo SB >>>
910    yqsurf = 0.0  ; yalb = 0.0 ; yalb_vis = 0.0
911!albedo SB <<<
912    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0   
913    ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yu1 = 0.0   
914    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
915    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
916    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0
917    yrugoro = 0.0 ; ywindsp = 0.0   
918!!    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
919    yfluxlat=0.0
920!!    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
921!!    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0
922    yqsol = 0.0   
923    ytherm = 0.0  ; ytke=0.
924    ! Martin
925    ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
926    yalb3_new = 0.0  ; ysissnow = 0.0
927    ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
928    ! Martin
929
930!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
931    ytke_x=0.     ; ytke_w=0.        ; ywake_dltke=0.
932    y_d_t_x=0.    ; y_d_t_w=0.       ; y_d_q_x=0.      ; y_d_q_w=0.
933!!    d_t_w=0.      ; d_q_w=0.         
934!!    d_t_x=0.      ; d_q_x=0.
935!!    d_wake_dlt=0.    ; d_wake_dlq=0.
936    yfluxlat_x=0. ; yfluxlat_w=0.
937    ywake_s=0.    ; ywake_cstar=0.   ;ywake_dens=0.
938!!!
939!!! nrlmd le 13/06/2011
940    tau_eq=0.     ; delta_coef=0.
941    y_delta_flux_t1=0.
942    ydtsurf_th=0.
943    yts_x=0.      ; yts_w=0.
944    y_delta_tsurf=0.
945!!!
946    ytsoil = 999999.
947
948
949! 2c) Initialization of all local variables computed within the subsurface loop and used later on
950!****************************************************************************************
951    d_t_diss_x(:,:) = 0. ;        d_t_diss_w(:,:) = 0.
952    d_u_x(:,:)=0. ;               d_u_w(:,:)=0.
953    d_v_x(:,:)=0. ;               d_v_w(:,:)=0.
954    flux_t_x(:,:,:)=0. ;          flux_t_w(:,:,:)=0.
955    flux_q_x(:,:,:)=0. ;          flux_q_w(:,:,:)=0.
956!
957!jyg<
958    flux_u_x(:,:,:)=0. ;          flux_u_w(:,:,:)=0.
959    flux_v_x(:,:,:)=0. ;          flux_v_w(:,:,:)=0.
960    fluxlat_x(:,:)=0. ;           fluxlat_w(:,:)=0.
961!>jyg
962!
963!jyg<
964! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces
965! actually present in the grid cell  ==> value set to 999999.
966!                           
967!jyg<
968       ustar(:,:)   = 999999.
969       wstar(:,:)   = 999999.
970       windsp(:,:)  = SQRT(u10m(:,:)**2 + v10m(:,:)**2 )
971       u10m(:,:)    = 999999.
972       v10m(:,:)    = 999999.
973!>jyg
974!
975       pblh(:,:)   = 999999.        ! Hauteur de couche limite
976       plcl(:,:)   = 999999.        ! Niveau de condensation de la CLA
977       capCL(:,:)  = 999999.        ! CAPE de couche limite
978       oliqCL(:,:) = 999999.        ! eau_liqu integree de couche limite
979       cteiCL(:,:) = 999999.        ! cloud top instab. crit. couche limite
980       pblt(:,:)   = 999999.        ! T a la Hauteur de couche limite
981       therm(:,:)  = 999999.
982       trmb1(:,:)  = 999999.        ! deep_cape
983       trmb2(:,:)  = 999999.        ! inhibition
984       trmb3(:,:)  = 999999.        ! Point Omega
985!
986       t2m_x(:,:)    = 999999.
987       q2m_x(:,:)    = 999999.
988       ustar_x(:,:)   = 999999.
989       wstar_x(:,:)   = 999999.
990       u10m_x(:,:)   = 999999.
991       v10m_x(:,:)   = 999999.
992!                           
993       pblh_x(:,:)   = 999999.      ! Hauteur de couche limite
994       plcl_x(:,:)   = 999999.      ! Niveau de condensation de la CLA
995       capCL_x(:,:)  = 999999.      ! CAPE de couche limite
996       oliqCL_x(:,:) = 999999.      ! eau_liqu integree de couche limite
997       cteiCL_x(:,:) = 999999.      ! cloud top instab. crit. couche limite
998       pblt_x(:,:)   = 999999.      ! T a la Hauteur de couche limite
999       therm_x(:,:)  = 999999.     
1000       trmb1_x(:,:)  = 999999.      ! deep_cape
1001       trmb2_x(:,:)  = 999999.      ! inhibition
1002       trmb3_x(:,:)  = 999999.      ! Point Omega
1003!
1004       t2m_w(:,:)    = 999999.
1005       q2m_w(:,:)    = 999999.
1006       ustar_w(:,:)   = 999999.
1007       wstar_w(:,:)   = 999999.
1008       u10m_w(:,:)   = 999999.
1009       v10m_w(:,:)   = 999999.
1010                           
1011       pblh_w(:,:)   = 999999.      ! Hauteur de couche limite
1012       plcl_w(:,:)   = 999999.      ! Niveau de condensation de la CLA
1013       capCL_w(:,:)  = 999999.      ! CAPE de couche limite
1014       oliqCL_w(:,:) = 999999.      ! eau_liqu integree de couche limite
1015       cteiCL_w(:,:) = 999999.      ! cloud top instab. crit. couche limite
1016       pblt_w(:,:)   = 999999.      ! T a la Hauteur de couche limite
1017       therm_w(:,:)  = 999999.     
1018       trmb1_w(:,:)  = 999999.      ! deep_cape
1019       trmb2_w(:,:)  = 999999.      ! inhibition
1020       trmb3_w(:,:)  = 999999.      ! Point Omega
1021!!!     
1022!
1023!!!
1024!****************************************************************************************
1025! 3) - Calculate pressure thickness of each layer
1026!    - Calculate the wind at first layer
1027!    - Mean calculations of albedo
1028!    - Calculate net radiance at sub-surface
1029!****************************************************************************************
1030    DO k = 1, klev
1031       DO i = 1, klon
1032          delp(i,k) = paprs(i,k)-paprs(i,k+1)
1033       ENDDO
1034    ENDDO
1035
1036!****************************************************************************************
1037! Test for rugos........ from physiq.. A la fin plutot???
1038!
1039!****************************************************************************************
1040
1041    DO nsrf = 1, nbsrf
1042       DO i = 1, klon
1043          z0m(i,nsrf) = MAX(z0m(i,nsrf),z0min)
1044          z0h(i,nsrf) = MAX(z0h(i,nsrf),z0min)
1045       ENDDO
1046    ENDDO
1047
1048! Mean calculations of albedo
1049!
1050! * alb  : mean albedo for whole SW interval
1051!
1052! Mean albedo for grid point
1053! * alb_m  : mean albedo at whole SW interval
1054
1055    alb_dir_m(:,:) = 0.0
1056    alb_dif_m(:,:) = 0.0
1057    DO k = 1, nsw
1058     DO nsrf = 1, nbsrf
1059       DO i = 1, klon
1060          alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf)
1061          alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf)
1062       ENDDO
1063     ENDDO
1064    ENDDO
1065
1066! We here suppose the fraction f1 of incoming radiance of visible radiance
1067! as a fraction of all shortwave radiance
1068    f1 = 0.5
1069!    f1 = 1    ! put f1=1 to recreate old calculations
1070
1071!f1 is already included with SFRWL values in each surf files
1072    alb=0.0
1073    DO k=1,nsw
1074      DO nsrf = 1, nbsrf
1075        DO i = 1, klon
1076            alb(i,nsrf) = alb(i,nsrf) + alb_dir(i,k,nsrf)*SFRWL(k)
1077        ENDDO
1078      ENDDO
1079    ENDDO
1080
1081    alb_m=0.0
1082    DO k = 1,nsw
1083      DO i = 1, klon
1084        alb_m(i) = alb_m(i) + alb_dir_m(i,k)*SFRWL(k)
1085      END DO
1086    ENDDO
1087!albedo SB <<<
1088
1089
1090
1091! Calculation of mean temperature at surface grid points
1092    ztsol(:) = 0.0
1093    DO nsrf = 1, nbsrf
1094       DO i = 1, klon
1095          ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
1096       ENDDO
1097    ENDDO
1098
1099! Linear distrubution on sub-surface of long- and shortwave net radiance
1100    DO nsrf = 1, nbsrf
1101       DO i = 1, klon
1102          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
1103
1104!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1105!         ! Martin
1106! Apparently introduced for sisvat but not used
1107!         sollwd(i,nsrf)= sollwd_m(i)
1108!         ! Martin
1109!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1110
1111          solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
1112       ENDDO
1113    ENDDO
1114!
1115!<al1: second order corrections
1116!- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
1117   IF (iflag_order2_sollw == 1) THEN
1118    meansqT(:) = 0. ! as working buffer
1119    DO nsrf = 1, nbsrf
1120     DO i = 1, klon
1121      meansqT(i) = meansqT(i)+(ts(i,nsrf)-ztsol(i))**2 *pctsrf(i,nsrf)
1122     END DO
1123    END DO
1124    DO nsrf = 1, nbsrf
1125     DO i = 1, klon
1126      sollw(i,nsrf) = sollw(i,nsrf) &
1127                + 6.0*RSIGMA*ztsol(i)**2 *(meansqT(i)-(ztsol(i)-ts(i,nsrf))**2)
1128     ENDDO
1129    ENDDO
1130   ENDIF   ! iflag_order2_sollw == 1
1131!>al1
1132
1133!****************************************************************************************
1134! 4) Loop over different surfaces
1135!
1136! Only points containing a fraction of the sub surface will be treated.
1137!
1138!****************************************************************************************
1139   
1140    loop_nbsrf: DO nsrf = 1, nbsrf
1141       IF (prt_level >=10) print *,' Loop nsrf ',nsrf
1142
1143! Search for index(ni) and size(knon) of domaine to treat
1144       ni(:) = 0
1145       knon  = 0
1146       DO i = 1, klon
1147          IF (pctsrf(i,nsrf) > 0.) THEN
1148             knon = knon + 1
1149             ni(knon) = i
1150          ENDIF
1151       ENDDO
1152
1153!!! jyg le 19/08/2012
1154!       IF (knon <= 0) THEN
1155!         IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf
1156!         cycle loop_nbsrf
1157!       ENDIF
1158!!!
1159
1160       ! write index, with IOIPSL
1161       IF (debugindex .AND. mpi_size==1) THEN
1162          tabindx(:)=0.
1163          DO i=1,knon
1164             tabindx(i)=REAL(i)
1165          END DO
1166          debugtab(:,:) = 0.
1167          ndexbg(:) = 0
1168          CALL gath2cpl(tabindx,debugtab,knon,ni)
1169          CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,nbp_lon*nbp_lat, ndexbg)
1170       ENDIF
1171       
1172!****************************************************************************************
1173! 5) Compress variables
1174!
1175!****************************************************************************************
1176
1177       DO j = 1, knon
1178          i = ni(j)
1179          ypct(j)    = pctsrf(i,nsrf)
1180          yts(j)     = ts(i,nsrf)
1181          ysnow(j)   = snow(i,nsrf)
1182          yqsurf(j)  = qsurf(i,nsrf)
1183          yalb(j)    = alb(i,nsrf)
1184!albedo SB >>>
1185          yalb_vis(j) = alb_dir(i,1,nsrf)
1186          if(nsw==6)then
1187            yalb_vis(j)=(alb_dir(i,1,nsrf)*SFRWL(1)+alb_dir(i,2,nsrf)*SFRWL(2) &
1188              +alb_dir(i,3,nsrf)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))
1189          endif
1190!albedo SB <<<
1191          yrain_f(j) = rain_f(i)
1192          ysnow_f(j) = snow_f(i)
1193          yagesno(j) = agesno(i,nsrf)
1194          yfder(j)   = fder(i)
1195          ylwdown(j) = lwdown_m(i)
1196          ygustiness(j) = gustiness(i)
1197          ysolsw(j)  = solsw(i,nsrf)
1198          ysollw(j)  = sollw(i,nsrf)
1199          yz0m(j)  = z0m(i,nsrf)
1200          yz0h(j)  = z0h(i,nsrf)
1201          yrugoro(j) = rugoro(i)
1202          yu1(j)     = u(i,1)
1203          yv1(j)     = v(i,1)
1204          ypaprs(j,klev+1) = paprs(i,klev+1)
1205!jyg<
1206!!          ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )
1207          ywindsp(j) = windsp(i,nsrf)
1208!>jyg
1209          ! Martin
1210          yzsig(j)   = zsig(i)
1211          ycldt(j)   = cldt(i)
1212          yrmu0(j)   = rmu0(i)
1213          ! Martin
1214!!! nrlmd le 13/06/2011
1215          y_delta_tsurf(j)=delta_tsurf(i,nsrf)
1216!!!
1217       END DO
1218
1219       DO k = 1, klev
1220          DO j = 1, knon
1221             i = ni(j)
1222             ypaprs(j,k) = paprs(i,k)
1223             ypplay(j,k) = pplay(i,k)
1224             ydelp(j,k)  = delp(i,k)
1225          ENDDO
1226       ENDDO
1227!!! jyg le 07/02/2012 et le 10/04/2013
1228        DO k = 1, klev
1229          DO j = 1, knon
1230             i = ni(j)
1231!jyg<
1232!!             ytke(j,k)   = tke(i,k,nsrf)
1233             ytke(j,k)   = tke_x(i,k,nsrf)
1234!>jyg
1235             yu(j,k) = u(i,k)
1236             yv(j,k) = v(i,k)
1237             yt(j,k) = t(i,k)
1238             yq(j,k) = q(i,k)
1239          ENDDO
1240        ENDDO
1241!
1242       IF (iflag_split .ge.1) THEN
1243!!! nrlmd le 02/05/2011
1244        DO k = 1, klev
1245          DO j = 1, knon
1246             i = ni(j)
1247             yu_x(j,k) = u(i,k)
1248             yv_x(j,k) = v(i,k)
1249             yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k)
1250             yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k)
1251             yu_w(j,k) = u(i,k)
1252             yv_w(j,k) = v(i,k)
1253             yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k)
1254             yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
1255!!!
1256          ENDDO
1257        ENDDO
1258!!! nrlmd le 02/05/2011
1259        DO k = 1, klev+1
1260          DO j = 1, knon
1261             i = ni(j)
1262!jyg<
1263!!             ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)
1264!!             ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)
1265!!             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
1266!!             ytke(j,k)     = tke(i,k,nsrf)
1267!
1268             ytke_x(j,k)      = tke_x(i,k,nsrf)
1269             ytke(j,k)        = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)
1270             ytke_w(j,k)      = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
1271             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
1272!>jyg
1273          ENDDO
1274        ENDDO
1275!!!
1276!!! jyg le 07/02/2012
1277        DO j = 1, knon
1278          i = ni(j)
1279          ywake_s(j)=wake_s(i)
1280          ywake_cstar(j)=wake_cstar(i)
1281          ywake_dens(j)=wake_dens(i)
1282        ENDDO
1283!!!
1284!!! nrlmd le 13/06/2011
1285        DO j=1,knon
1286         yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j)
1287         yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j)
1288        ENDDO
1289!!!
1290       ENDIF  ! (iflag_split .ge.1)
1291!!!
1292       DO k = 1, nsoilmx
1293          DO j = 1, knon
1294             i = ni(j)
1295             ytsoil(j,k) = ftsoil(i,k,nsrf)
1296          END DO
1297       END DO
1298       
1299       ! qsol(water height in soil) only for bucket continental model
1300       IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
1301          DO j = 1, knon
1302             i = ni(j)
1303             yqsol(j) = qsol(i)
1304          END DO
1305       ENDIF
1306       
1307!****************************************************************************************
1308! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
1309!
1310!****************************************************************************************
1311
1312!!! jyg le 07/02/2012
1313       IF (iflag_split .eq.0) THEN
1314!!!
1315!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1316! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1317!       CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1318!           yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
1319!           yts, yqsurf, yrugos, &
1320!           ycdragm, ycdragh )
1321! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1322        DO i = 1, knon
1323!          print*,'PBL ',i,RD
1324!          print*,'PBL ',yt(i,1),ypaprs(i,1),ypplay(i,1)
1325           zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1326                * (ypaprs(i,1)-ypplay(i,1))
1327           speed(i) = SQRT(yu(i,1)**2+yv(i,1)**2)
1328        END DO
1329        CALL cdrag(knon, nsrf, &
1330            speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1),&
1331            yts, yqsurf, yz0m, yz0h, &
1332            ycdragm, ycdragh, zri1, pref )
1333
1334! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
1335     IF (ok_prescr_ust) then
1336      DO i = 1, knon
1337       print *,'ycdragm avant=',ycdragm(i)
1338       vent= sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))
1339!      ycdragm(i) = ust*ust/(1.+(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)))
1340!      ycdragm(i) = ust*ust/((1.+sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))) &
1341!     *sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)))
1342       ycdragm(i) = ust*ust/(1.+vent)/vent
1343!      print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1)
1344      ENDDO
1345     ENDIF
1346        IF (prt_level >=10) print *,'clcdrag -> ycdragh ', ycdragh
1347       ELSE  !(iflag_split .eq.0)
1348
1349! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1350!       CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1351!           yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &
1352!           yts_x, yqsurf, yrugos, &
1353!           ycdragm_x, ycdragh_x )
1354! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1355        DO i = 1, knon
1356           zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1357                * (ypaprs(i,1)-ypplay(i,1))
1358           speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2)
1359        END DO
1360        CALL cdrag(knon, nsrf, &
1361            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),&
1362            yts_x, yqsurf, yz0m, yz0h, &
1363            ycdragm_x, ycdragh_x, zri1_x, pref_x )
1364
1365! --- special Dice. JYG+MPL 25112013
1366        IF (ok_prescr_ust) then
1367         DO i = 1, knon
1368!         print *,'ycdragm_x avant=',ycdragm_x(i)
1369          vent= sqrt(yu_x(i,1)*yu_x(i,1)+yv_x(i,1)*yv_x(i,1))
1370          ycdragm_x(i) = ust*ust/(1.+vent)/vent
1371!         print *,'ycdragm_x ust yu yv apres=',ycdragm_x(i),ust,yu_x(i,1),yv_x(i,1)
1372         ENDDO
1373        ENDIF
1374        IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x
1375!
1376! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1377!        CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1378!            yu_w(:,1), yv_w(:,1), yt_w(:,1), yq_w(:,1), &
1379!            yts_w, yqsurf, yz0m, &
1380!            ycdragm_w, ycdragh_w )
1381! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1382        DO i = 1, knon
1383           zgeo1_w(i) = RD * yt_w(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1384                * (ypaprs(i,1)-ypplay(i,1))
1385           speed_w(i) = SQRT(yu_w(i,1)**2+yv_w(i,1)**2)
1386        END DO
1387        CALL cdrag(knon, nsrf, &
1388            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),&
1389            yts_w, yqsurf, yz0m, yz0h, &
1390            ycdragm_w, ycdragh_w, zri1_w, pref_w )
1391
1392! --- special Dice. JYG+MPL 25112013 puis BOMEX
1393        IF (ok_prescr_ust) then
1394         DO i = 1, knon
1395!         print *,'ycdragm_w avant=',ycdragm_w(i)
1396          vent= sqrt(yu_w(i,1)*yu_w(i,1)+yv_w(i,1)*yv_w(i,1))
1397          ycdragm_w(i) = ust*ust/(1.+vent)/vent
1398!         print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1)
1399         ENDDO
1400        ENDIF
1401        IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w
1402!!!
1403       ENDIF  ! (iflag_split .eq.0)
1404!!!
1405       
1406
1407!****************************************************************************************
1408! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.
1409!
1410!****************************************************************************************
1411
1412!!! jyg le 07/02/2012
1413       IF (iflag_split .eq.0) THEN
1414!!!
1415!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1416      IF (prt_level >=10) THEN
1417      print *,' args coef_diff_turb: yu ',  yu 
1418      print *,' args coef_diff_turb: yv ',  yv 
1419      print *,' args coef_diff_turb: yq ',  yq 
1420      print *,' args coef_diff_turb: yt ',  yt 
1421      print *,' args coef_diff_turb: yts ', yts 
1422      print *,' args coef_diff_turb: yz0m ', yz0m 
1423      print *,' args coef_diff_turb: yqsurf ', yqsurf 
1424      print *,' args coef_diff_turb: ycdragm ', ycdragm
1425      print *,' args coef_diff_turb: ycdragh ', ycdragh
1426      print *,' args coef_diff_turb: ytke ', ytke
1427       ENDIF
1428        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
1429            ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
1430            ycoefm, ycoefh, ytke)
1431       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1432! In this case, coef_diff_turb is called for the Cd only
1433       DO k = 2, klev
1434          DO j = 1, knon
1435             i = ni(j)
1436             ycoefh(j,k)   = zcoefh(i,k,nsrf)
1437             ycoefm(j,k)   = zcoefm(i,k,nsrf)
1438          ENDDO
1439       ENDDO
1440       ENDIF
1441        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh
1442!
1443       ELSE  !(iflag_split .eq.0)
1444      IF (prt_level >=10) THEN
1445      print *,' args coef_diff_turb: yu_x ',  yu_x 
1446      print *,' args coef_diff_turb: yv_x ',  yv_x 
1447      print *,' args coef_diff_turb: yq_x ',  yq_x 
1448      print *,' args coef_diff_turb: yt_x ',  yt_x 
1449      print *,' args coef_diff_turb: yts_x ', yts_x 
1450      print *,' args coef_diff_turb: yqsurf ', yqsurf 
1451      print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x
1452      print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x
1453      print *,' args coef_diff_turb: ytke_x ', ytke_x
1454       ENDIF
1455        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
1456            ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf, ycdragm_x, &
1457            ycoefm_x, ycoefh_x, ytke_x)
1458       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1459! In this case, coef_diff_turb is called for the Cd only
1460       DO k = 2, klev
1461          DO j = 1, knon
1462             i = ni(j)
1463             ycoefh_x(j,k)   = zcoefh(i,k,nsrf)
1464             ycoefm_x(j,k)   = zcoefm(i,k,nsrf)
1465          ENDDO
1466       ENDDO
1467       ENDIF
1468        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x
1469!
1470      IF (prt_level >=10) THEN
1471      print *,' args coef_diff_turb: yu_w ',  yu_w 
1472      print *,' args coef_diff_turb: yv_w ',  yv_w 
1473      print *,' args coef_diff_turb: yq_w ',  yq_w 
1474      print *,' args coef_diff_turb: yt_w ',  yt_w 
1475      print *,' args coef_diff_turb: yts_w ', yts_w 
1476      print *,' args coef_diff_turb: yqsurf ', yqsurf 
1477      print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w
1478      print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w
1479      print *,' args coef_diff_turb: ytke_w ', ytke_w
1480       ENDIF
1481        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
1482            ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf, ycdragm_w, &
1483            ycoefm_w, ycoefh_w, ytke_w)
1484       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
1485! In this case, coef_diff_turb is called for the Cd only
1486       DO k = 2, klev
1487          DO j = 1, knon
1488             i = ni(j)
1489             ycoefh_w(j,k)   = zcoefh(i,k,nsrf)
1490             ycoefm_w(j,k)   = zcoefm(i,k,nsrf)
1491          ENDDO
1492       ENDDO
1493       ENDIF
1494        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w
1495!
1496!!!jyg le 10/04/2013
1497!!   En attendant de traiter le transport des traceurs dans les poches froides, formule
1498!!   arbitraire pour ycoefh et ycoefm
1499      DO k = 2,klev
1500        DO j = 1,knon
1501         ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))
1502         ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))
1503        ENDDO
1504      ENDDO
1505!!!
1506       ENDIF  ! (iflag_split .eq.0)
1507!!!
1508       
1509!****************************************************************************************
1510!
1511! 8) "La descente" - "The downhill"
1512
1513!  climb_hq_down and climb_wind_down calculate the coefficients
1514!  Ccoef_X et Dcoef_X for X=[H, Q, U, V].
1515!  Only the coefficients at surface for H and Q are returned.
1516!
1517!****************************************************************************************
1518
1519! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
1520!!! jyg le 07/02/2012
1521       IF (iflag_split .eq.0) THEN
1522!!!
1523!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1524        CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
1525            ydelp, yt, yq, dtime, &
1526!!! jyg le 09/05/2011
1527            CcoefH, CcoefQ, DcoefH, DcoefQ, &
1528            Kcoef_hq, gama_q, gama_h, &
1529!!!
1530            AcoefH, AcoefQ, BcoefH, BcoefQ)
1531       ELSE  !(iflag_split .eq.0)
1532        CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
1533            ydelp, yt_x, yq_x, dtime, &
1534!!! nrlmd le 02/05/2011
1535            CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
1536            Kcoef_hq_x, gama_q_x, gama_h_x, &
1537!!!
1538            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x)
1539!
1540        CALL climb_hq_down(knon, ycoefh_w, ypaprs, ypplay, &
1541            ydelp, yt_w, yq_w, dtime, &
1542!!! nrlmd le 02/05/2011
1543            CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
1544            Kcoef_hq_w, gama_q_w, gama_h_w, &
1545!!!
1546            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w)
1547!!!
1548       ENDIF  ! (iflag_split .eq.0)
1549!!!
1550
1551! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
1552!!! jyg le 07/02/2012
1553       IF (iflag_split .eq.0) THEN
1554!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1555        CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
1556!!! jyg le 09/05/2011
1557            CcoefU, CcoefV, DcoefU, DcoefV, &
1558            Kcoef_m, alf_1, alf_2, &
1559!!!
1560            AcoefU, AcoefV, BcoefU, BcoefV)
1561       ELSE  ! (iflag_split .eq.0)
1562        CALL climb_wind_down(knon, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &
1563!!! nrlmd le 02/05/2011
1564            CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
1565            Kcoef_m_x, alf_1_x, alf_2_x, &
1566!!!
1567            AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)
1568!
1569        CALL climb_wind_down(knon, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &
1570!!! nrlmd le 02/05/2011
1571            CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
1572            Kcoef_m_w, alf_1_w, alf_2_w, &
1573!!!
1574            AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)
1575!!!     
1576       ENDIF  ! (iflag_split .eq.0)
1577!!!
1578
1579!****************************************************************************************
1580! 9) Small calculations
1581!
1582!****************************************************************************************
1583
1584! - Reference pressure is given the values at surface level         
1585       ypsref(:) = ypaprs(:,1) 
1586
1587! - CO2 field on 2D grid to be sent to ORCHIDEE
1588!   Transform to compressed field
1589       IF (carbon_cycle_cpl) THEN
1590          DO i=1,knon
1591             r_co2_ppm(i) = co2_send(ni(i))
1592          END DO
1593       ELSE
1594          r_co2_ppm(:) = co2_ppm     ! Constant field
1595       END IF
1596
1597!!! nrlmd le 13/06/2011
1598!----- On finit le calcul des coefficients d'\'echange:on multiplie le cdrag par le module du vent et la densit\'e dans la premi\`ere couche
1599!          Kech_h_x(j) = ycdragh_x(j) * &
1600!             (1.0+SQRT(yu_x(j,1)**2+yv_x(j,1)**2)) * &
1601!             ypplay(j,1)/(RD*yt_x(j,1))
1602!          Kech_h_w(j) = ycdragh_w(j) * &
1603!             (1.0+SQRT(yu_w(j,1)**2+yv_w(j,1)**2)) * &
1604!             ypplay(j,1)/(RD*yt_w(j,1))
1605!          Kech_h(j) = (1.-ywake_s(j))*Kech_h_x(j)+ywake_s(j)*Kech_h_w(j)
1606!
1607!          Kech_m_x(j) = ycdragm_x(j) * &
1608!             (1.0+SQRT(yu_x(j,1)**2+yv_x(j,1)**2)) * &
1609!             ypplay(j,1)/(RD*yt_x(j,1))
1610!          Kech_m_w(j) = ycdragm_w(j) * &
1611!             (1.0+SQRT(yu_w(j,1)**2+yv_w(j,1)**2)) * &
1612!             ypplay(j,1)/(RD*yt_w(j,1))
1613!          Kech_m(j) = (1.-ywake_s(j))*Kech_m_x(j)+ywake_s(j)*Kech_m_w(j)
1614!!!
1615
1616!!! nrlmd le 02/05/2011  -----------------------On raccorde les 2 colonnes dans la couche 1
1617!----------------------------------------------------------------------------------------
1618!!! jyg le 07/02/2012
1619       IF (iflag_split .eq.1) THEN
1620!!!
1621!!! jyg le 09/04/2013 ; passage aux nouvelles expressions en differences
1622
1623        DO j=1,knon
1624!
1625! Calcul des coefficients d echange
1626         mod_wind_x = 1.0+SQRT(yu_x(j,1)**2+yv_x(j,1)**2)
1627         mod_wind_w = 1.0+SQRT(yu_w(j,1)**2+yv_w(j,1)**2)
1628         rho1 = ypplay(j,1)/(RD*yt(j,1))
1629         Kech_h_x(j) = ycdragh_x(j) * mod_wind_x * rho1
1630         Kech_h_w(j) = ycdragh_w(j) * mod_wind_w * rho1
1631         Kech_m_x(j) = ycdragm_x(j) * mod_wind_x * rho1
1632         Kech_m_w(j) = ycdragm_w(j) * mod_wind_w * rho1
1633!
1634         dd_Kh = Kech_h_w(j) - Kech_h_x(j)
1635         dd_Km = Kech_m_w(j) - Kech_m_x(j)
1636         IF (prt_level >=10) THEN
1637          print *,' mod_wind_x, mod_wind_w ', mod_wind_x, mod_wind_w
1638          print *,' rho1 ',rho1
1639          print *,' ycdragh_x(j),ycdragm_x(j) ',ycdragh_x(j),ycdragm_x(j)
1640          print *,' ycdragh_w(j),ycdragm_w(j) ',ycdragh_w(j),ycdragm_w(j)
1641          print *,' dd_Kh: ',dd_KH
1642         ENDIF
1643!
1644         Kech_h(j) = Kech_h_x(j) + ywake_s(j)*dd_Kh
1645         Kech_m(j) = Kech_m_x(j) + ywake_s(j)*dd_Km
1646!
1647! Calcul des coefficients d echange corriges des retroactions
1648        Kech_H_xp(j) = Kech_h_x(j)/(1.-BcoefH_x(j)*Kech_h_x(j)*dtime)
1649        Kech_H_wp(j) = Kech_h_w(j)/(1.-BcoefH_w(j)*Kech_h_w(j)*dtime)
1650        Kech_Q_xp(j) = Kech_h_x(j)/(1.-BcoefQ_x(j)*Kech_h_x(j)*dtime)
1651        Kech_Q_wp(j) = Kech_h_w(j)/(1.-BcoefQ_w(j)*Kech_h_w(j)*dtime)
1652        Kech_U_xp(j) = Kech_m_x(j)/(1.-BcoefU_x(j)*Kech_m_x(j)*dtime)
1653        Kech_U_wp(j) = Kech_m_w(j)/(1.-BcoefU_w(j)*Kech_m_w(j)*dtime)
1654        Kech_V_xp(j) = Kech_m_x(j)/(1.-BcoefV_x(j)*Kech_m_x(j)*dtime)
1655        Kech_V_wp(j) = Kech_m_w(j)/(1.-BcoefV_w(j)*Kech_m_w(j)*dtime)
1656!
1657         dd_KHp = Kech_H_wp(j) - Kech_H_xp(j)
1658         dd_KQp = Kech_Q_wp(j) - Kech_Q_xp(j)
1659         dd_KUp = Kech_U_wp(j) - Kech_U_xp(j)
1660         dd_KVp = Kech_V_wp(j) - Kech_V_xp(j)
1661!
1662        Kech_Hp(j) = Kech_H_xp(j) + ywake_s(j)*dd_KHp
1663        Kech_Qp(j) = Kech_Q_xp(j) + ywake_s(j)*dd_KQp
1664        Kech_Up(j) = Kech_U_xp(j) + ywake_s(j)*dd_KUp
1665        Kech_Vp(j) = Kech_V_xp(j) + ywake_s(j)*dd_KVp
1666!
1667! Calcul des differences w-x
1668       dd_CM = ycdragm_w(j) - ycdragm_x(j)
1669       dd_CH = ycdragh_w(j) - ycdragh_x(j)
1670       dd_u = yu_w(j,1) - yu_x(j,1)
1671       dd_v = yv_w(j,1) - yv_x(j,1)
1672       dd_t = yt_w(j,1) - yt_x(j,1)
1673       dd_q = yq_w(j,1) - yq_x(j,1)
1674       dd_AH = AcoefH_w(j) - AcoefH_x(j)
1675       dd_AQ = AcoefQ_w(j) - AcoefQ_x(j)
1676       dd_AU = AcoefU_w(j) - AcoefU_x(j)
1677       dd_AV = AcoefV_w(j) - AcoefV_x(j)
1678       dd_BH = BcoefH_w(j) - BcoefH_x(j)
1679       dd_BQ = BcoefQ_w(j) - BcoefQ_x(j)
1680       dd_BU = BcoefU_w(j) - BcoefU_x(j)
1681       dd_BV = BcoefV_w(j) - BcoefV_x(j)
1682!
1683       IF (prt_level >=10) THEN
1684          print *,'Variables pour la fusion : Kech_H_xp(j)' ,Kech_H_xp(j)
1685          print *,'Variables pour la fusion : Kech_H_wp(j)' ,Kech_H_wp(j)
1686          print *,'Variables pour la fusion : Kech_Hp(j)' ,Kech_Hp(j)
1687          print *,'Variables pour la fusion : Kech_h(j)' ,Kech_h(j)
1688       ENDIF
1689!
1690! Calcul des coef A, B \'equivalents dans la couche 1
1691!
1692       AcoefH(j) = AcoefH_x(j) + ywake_s(j)*(Kech_H_wp(j)/Kech_Hp(j))*dd_AH
1693       AcoefQ(j) = AcoefQ_x(j) + ywake_s(j)*(Kech_Q_wp(j)/Kech_Qp(j))*dd_AQ
1694       AcoefU(j) = AcoefU_x(j) + ywake_s(j)*(Kech_U_wp(j)/Kech_Up(j))*dd_AU
1695       AcoefV(j) = AcoefV_x(j) + ywake_s(j)*(Kech_V_wp(j)/Kech_Vp(j))*dd_AV
1696!
1697       BcoefH(j) = BcoefH_x(j) + ywake_s(j)*BcoefH_x(j)*(dd_Kh/Kech_h(j))*(1.+Kech_H_wp(j)/Kech_Hp(j)) &
1698                               + ywake_s(j)*(Kech_H_wp(j)/Kech_Hp(j))*(Kech_h_w(j)/Kech_h(j))*dd_BH
1699
1700       BcoefQ(j) = BcoefQ_x(j) + ywake_s(j)*BcoefQ_x(j)*(dd_Kh/Kech_h(j))*(1.+Kech_Q_wp(j)/Kech_Qp(j)) &
1701                               + ywake_s(j)*(Kech_Q_wp(j)/Kech_Qp(j))*(Kech_h_w(j)/Kech_h(j))*dd_BQ
1702
1703       BcoefU(j) = BcoefU_x(j) + ywake_s(j)*BcoefU_x(j)*(dd_Km/Kech_h(j))*(1.+Kech_U_wp(j)/Kech_Up(j)) &
1704                               + ywake_s(j)*(Kech_U_wp(j)/Kech_Up(j))*(Kech_m_w(j)/Kech_m(j))*dd_BU
1705
1706       BcoefV(j) = BcoefV_x(j) + ywake_s(j)*BcoefV_x(j)*(dd_Km/Kech_h(j))*(1.+Kech_V_wp(j)/Kech_Vp(j)) &
1707                               + ywake_s(j)*(Kech_V_wp(j)/Kech_Vp(j))*(Kech_m_w(j)/Kech_m(j))*dd_BV
1708
1709!
1710! Calcul des cdrag \'equivalents dans la couche
1711!
1712       ycdragm(j) = ycdragm_x(j) + ywake_s(j)*dd_CM
1713       ycdragh(j) = ycdragh_x(j) + ywake_s(j)*dd_CH
1714!
1715! Calcul de T, q, u et v \'equivalents dans la couche 1
1716       yt(j,1) = yt_x(j,1) + ywake_s(j)*(Kech_h_w(j)/Kech_h(j))*dd_t
1717       yq(j,1) = yq_x(j,1) + ywake_s(j)*(Kech_h_w(j)/Kech_h(j))*dd_q
1718       yu(j,1) = yu_x(j,1) + ywake_s(j)*(Kech_m_w(j)/Kech_m(j))*dd_u
1719       yv(j,1) = yv_x(j,1) + ywake_s(j)*(Kech_m_w(j)/Kech_m(j))*dd_v
1720
1721
1722        ENDDO
1723!!!
1724       ENDIF  ! (iflag_split .eq.1)
1725!!!
1726
1727!****************************************************************************************
1728!
1729! Calulate t2m and q2m for the case of calculation at land grid points
1730! t2m and q2m are needed as input to ORCHIDEE
1731!
1732!****************************************************************************************
1733       IF (nsrf == is_ter) THEN
1734
1735          DO i = 1, knon
1736             zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1737                  * (ypaprs(i,1)-ypplay(i,1))
1738          END DO
1739
1740          ! Calculate the temperature et relative humidity at 2m and the wind at 10m
1741          CALL stdlevvar(klon, knon, is_ter, zxli, &
1742               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
1743               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
1744               yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1745         
1746       END IF
1747
1748!****************************************************************************************
1749!
1750! 10) Switch according to current surface
1751!     It is necessary to start with the continental surfaces because the ocean
1752!     needs their run-off.
1753!
1754!****************************************************************************************
1755       SELECT CASE(nsrf)
1756     
1757       CASE(is_ter)
1758!          print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
1759          CALL surf_land(itap, dtime, date0, jour, knon, ni,&
1760               rlon, rlat, yrmu0, &
1761               debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
1762               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1763               AcoefH, AcoefQ, BcoefH, BcoefQ, &
1764               AcoefU, AcoefV, BcoefU, BcoefV, &
1765               ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1766               ylwdown, yq2m, yt2m, &
1767               ysnow, yqsol, yagesno, ytsoil, &
1768               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1769               yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
1770               y_flux_u1, y_flux_v1 )
1771               
1772! Special DICE MPL 05082013 puis BOMEX
1773       IF (ok_prescr_ust) THEN
1774          do j=1,knon
1775!         ysnow(:)=0.
1776!         yqsol(:)=0.
1777!         yagesno(:)=50.
1778!         ytsoil(:,:)=300.
1779!         yz0_new(:)=0.001
1780!         yevap(:)=flat/RLVTT
1781!         yfluxlat(:)=-flat
1782!         yfluxsens(:)=-fsens
1783!         yqsurf(:)=0.
1784!         ytsurf_new(:)=tg
1785!         y_dflux_t(:)=0.
1786!         y_dflux_q(:)=0.
1787          y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
1788          y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
1789          enddo
1790      ENDIF
1791
1792     
1793       CASE(is_lic)
1794          ! Martin
1795          CALL surf_landice(itap, dtime, knon, ni, &
1796               rlon, rlat, debut, lafin, &
1797               yrmu0, ylwdown, yalb, ypphi(:,1), &
1798               ysolsw, ysollw, yts, ypplay(:,1), &
1799               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1800               AcoefH, AcoefQ, BcoefH, BcoefQ, &
1801               AcoefU, AcoefV, BcoefU, BcoefV, &
1802               ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1803               ysnow, yqsurf, yqsol, yagesno, &
1804               ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
1805               ytsurf_new, y_dflux_t, y_dflux_q, &
1806               yzsig, ycldt, &
1807               ysnowhgt, yqsnow, ytoice, ysissnow, &
1808               yalb3_new, yrunoff, &
1809               y_flux_u1, y_flux_v1)
1810
1811!jyg<
1812!!          alb3_lic(:)=0.
1813!>jyg
1814          DO j = 1, knon
1815             i = ni(j)
1816             alb3_lic(i) = yalb3_new(j)
1817             snowhgt(i)   = ysnowhgt(j)
1818             qsnow(i)     = yqsnow(j)
1819             to_ice(i)    = ytoice(j)
1820             sissnow(i)   = ysissnow(j)
1821             runoff(i)    = yrunoff(j)
1822          END DO
1823          ! Martin
1824! Special DICE MPL 05082013 puis BOMEX MPL 20150410
1825       IF (ok_prescr_ust) THEN
1826          DO j=1,knon
1827          y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
1828          y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
1829          ENDDO
1830      ENDIF
1831         
1832       CASE(is_oce)
1833           CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
1834               ywindsp, rmu0, yfder, yts, &
1835               itap, dtime, jour, knon, ni, &
1836               ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1837               AcoefH, AcoefQ, BcoefH, BcoefQ, &
1838               AcoefU, AcoefV, BcoefU, BcoefV, &
1839               ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
1840               ysnow, yqsurf, yagesno, &
1841               yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1842               ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
1843               y_flux_u1, y_flux_v1)
1844      IF (prt_level >=10) THEN
1845          print *,'arg de surf_ocean: ycdragh ',ycdragh
1846          print *,'arg de surf_ocean: ycdragm ',ycdragm
1847          print *,'arg de surf_ocean: yt ', yt
1848          print *,'arg de surf_ocean: yq ', yq
1849          print *,'arg de surf_ocean: yts ', yts
1850          print *,'arg de surf_ocean: AcoefH ',AcoefH
1851          print *,'arg de surf_ocean: AcoefQ ',AcoefQ
1852          print *,'arg de surf_ocean: BcoefH ',BcoefH
1853          print *,'arg de surf_ocean: BcoefQ ',BcoefQ
1854          print *,'arg de surf_ocean: yevap ',yevap
1855          print *,'arg de surf_ocean: yfluxsens ',yfluxsens
1856          print *,'arg de surf_ocean: yfluxlat ',yfluxlat
1857          print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new
1858       ENDIF
1859! Special DICE MPL 05082013 puis BOMEX MPL 20150410
1860       IF (ok_prescr_ust) THEN
1861          DO j=1,knon
1862          y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
1863          y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
1864          ENDDO
1865      ENDIF
1866         
1867       CASE(is_sic)
1868          CALL surf_seaice( &
1869!albedo SB >>>
1870               rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &
1871!albedo SB <<<
1872               itap, dtime, jour, knon, ni, &
1873               lafin, &
1874               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
1875               AcoefH, AcoefQ, BcoefH, BcoefQ, &
1876               AcoefU, AcoefV, BcoefU, BcoefV, &
1877               ypsref, yu1, yv1, ygustiness, pctsrf, &
1878               ysnow, yqsurf, yqsol, yagesno, ytsoil, &
1879!albedo SB >>>
1880               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
1881!albedo SB <<<
1882               ytsurf_new, y_dflux_t, y_dflux_q, &
1883               y_flux_u1, y_flux_v1)
1884         
1885! Special DICE MPL 05082013 puis BOMEX MPL 20150410
1886       IF (ok_prescr_ust) THEN
1887          DO j=1,knon
1888          y_flux_u1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yu(j,1)*ypplay(j,1)/RD/yt(j,1)
1889          y_flux_v1(j)=ycdragm(j)*(1.+sqrt(yu(j,1)*yu(j,1)+yv(j,1)*yv(j,1)))*yv(j,1)*ypplay(j,1)/RD/yt(j,1)
1890          ENDDO
1891      ENDIF
1892
1893       CASE DEFAULT
1894          WRITE(lunout,*) 'Surface index = ', nsrf
1895          abort_message = 'Surface index not valid'
1896          CALL abort_physic(modname,abort_message,1)
1897       END SELECT
1898
1899
1900!****************************************************************************************
1901! 11) - Calcul the increment of surface temperature
1902!
1903!****************************************************************************************
1904
1905       if (evap0>=0.) then
1906          yevap(:)=evap0
1907          yevap(:)=RLVTT*evap0
1908       endif
1909
1910
1911       y_d_ts(1:knon)   = ytsurf_new(1:knon) - yts(1:knon)
1912 
1913!****************************************************************************************
1914!
1915! 12) "La remontee" - "The uphill"
1916!
1917!  The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
1918!  for X=H, Q, U and V, for all vertical levels.
1919!
1920!****************************************************************************************
1921
1922!!!
1923!!! jyg le 10/04/2013
1924!!!
1925        IF (ok_flux_surf) THEN
1926          IF (prt_level >=10) THEN
1927           PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
1928          ENDIF
1929          y_flux_t1(:) =  fsens
1930          y_flux_q1(:) =  flat/RLVTT
1931          yfluxlat(:) =  flat
1932!
1933          IF (iflag_split .eq.0) THEN
1934             do j=1,knon
1935             Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)**2+yv(j,1)**2)) * &
1936                  ypplay(j,1)/(RD*yt(j,1))
1937             enddo
1938          ENDIF ! (iflag_split .eq.0)
1939
1940          DO j = 1, knon
1941            yt1_new=(1./RCPD)*(AcoefH(j)+BcoefH(j)*y_flux_t1(j)*dtime)
1942            ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)
1943          ENDDO
1944
1945          do j=1,knon
1946          y_d_ts(j) = ytsurf_new(j) - yts(j)
1947          enddo
1948
1949        ELSE ! (ok_flux_surf)
1950          do j=1,knon
1951          y_flux_t1(j) =  yfluxsens(j)
1952          y_flux_q1(j) = -yevap(j)
1953          enddo
1954        ENDIF
1955
1956       IF (prt_level >=10) THEN
1957        DO j=1,knon
1958         print*,'y_flux_t1,yfluxlat,wakes' &
1959 &             ,  y_flux_t1(j), yfluxlat(j), ywake_s(j)
1960         print*,'beta,ytsurf_new', ybeta(j), ytsurf_new(j)
1961         print*,'effusivity,facteur,cstar', effusivity, facteur,wake_cstar(j)
1962        ENDDO
1963       ENDIF
1964
1965!!! jyg le 07/02/2012 puis le 10/04/2013
1966       IF (iflag_split .eq.1) THEN
1967!!!
1968        DO j=1,knon
1969         y_delta_flux_t1(j) = ( Kech_H_wp(j)*Kech_H_xp(j)*(AcoefH_w(j)-AcoefH_x(j)) + &
1970                                y_flux_t1(j)*(Kech_H_wp(j)-Kech_H_xp(j)) ) / Kech_Hp(j)
1971         y_delta_flux_q1(j) = ( Kech_Q_wp(j)*Kech_Q_xp(j)*(AcoefQ_w(j)-AcoefQ_x(j)) + &
1972                                y_flux_q1(j)*(Kech_Q_wp(j)-Kech_Q_xp(j)) ) / Kech_Qp(j)
1973         y_delta_flux_u1(j) = ( Kech_U_wp(j)*Kech_U_xp(j)*(AcoefU_w(j)-AcoefU_x(j)) + &
1974                                y_flux_u1(j)*(Kech_U_wp(j)-Kech_U_xp(j)) ) / Kech_Up(j)
1975         y_delta_flux_v1(j) = ( Kech_V_wp(j)*Kech_V_xp(j)*(AcoefV_w(j)-AcoefV_x(j)) + &
1976                                y_flux_v1(j)*(Kech_V_wp(j)-Kech_V_xp(j)) ) / Kech_Vp(j)
1977!
1978         y_flux_t1_x(j)=y_flux_t1(j) - ywake_s(j)*y_delta_flux_t1(j)
1979         y_flux_t1_w(j)=y_flux_t1(j) + (1.-ywake_s(j))*y_delta_flux_t1(j)
1980         y_flux_q1_x(j)=y_flux_q1(j) - ywake_s(j)*y_delta_flux_q1(j)
1981         y_flux_q1_w(j)=y_flux_q1(j) + (1.-ywake_s(j))*y_delta_flux_q1(j)
1982         y_flux_u1_x(j)=y_flux_u1(j) - ywake_s(j)*y_delta_flux_u1(j)
1983         y_flux_u1_w(j)=y_flux_u1(j) + (1.-ywake_s(j))*y_delta_flux_u1(j)
1984         y_flux_v1_x(j)=y_flux_v1(j) - ywake_s(j)*y_delta_flux_v1(j)
1985         y_flux_v1_w(j)=y_flux_v1(j) + (1.-ywake_s(j))*y_delta_flux_v1(j)
1986!
1987         yfluxlat_x(j)=y_flux_q1_x(j)*RLVTT
1988         yfluxlat_w(j)=y_flux_q1_w(j)*RLVTT
1989
1990        ENDDO
1991!
1992 
1993!!jyg!!   A reprendre apres reflexion   ===============================================
1994!!jyg!!
1995!!jyg!!        DO j=1,knon
1996!!jyg!!!!! nrlmd le 13/06/2011
1997!!jyg!!
1998!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
1999!!jyg!!       IF (nsrf.eq.is_ter) THEN
2000!!jyg!!!----Calcul du coefficient delta_coeff
2001!!jyg!!          tau_eq(j)=(ywake_s(j)/2.)*(1./max(wake_cstar(j),0.01))*sqrt(0.4/(3.14*max(wake_dens(j),8e-12)))
2002!!jyg!!
2003!!jyg!!!          delta_coef(j)=dtime/(effusivity*sqrt(tau_eq(j)))
2004!!jyg!!          delta_coef(j)=facteur*sqrt(tau_eq(j))/effusivity
2005!!jyg!!!          delta_coef(j)=0.
2006!!jyg!!       ELSE
2007!!jyg!!         delta_coef(j)=0.
2008!!jyg!!       ENDIF
2009!!jyg!!
2010!!jyg!!!----Calcul de delta_tsurf
2011!!jyg!!         y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
2012!!jyg!!
2013!!jyg!!!----Si il n'y a pas des poches...
2014!!jyg!!         IF (wake_cstar(j).le.0.01) THEN
2015!!jyg!!           y_delta_tsurf(j)=0.
2016!!jyg!!           y_delta_flux_t1(j)=0.
2017!!jyg!!         ENDIF
2018!!jyg!!
2019!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)
2020!!jyg!!!!!!! jyg le 23/02/2012
2021!!jyg!!!!!!!
2022!!jyg!!!!        ybeta(j)=y_flux_q1(j)   /    &
2023!!jyg!!!! &        (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
2024!!jyg!!!!!!        ybeta(j)=-1.*yevap(j)   /    &
2025!!jyg!!!!!! &        (ywake_s(j)*Kech_h_w(j)*(yq_w(j,1)-yqsatsurf_w(j))+(1.-ywake_s(j))*Kech_h_x(j)*(yq_x(j,1)-yqsatsurf_x(j)))
2026!!jyg!!!!!!! fin jyg
2027!!jyg!!!!!
2028!!jyg!!
2029!!jyg!!       ENDDO
2030!!jyg!!
2031!!jyg!!!!! fin nrlmd le 13/06/2011
2032!!jyg!!
2033       IF (prt_level >=10) THEN
2034        DO j = 1, knon
2035         print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
2036         print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
2037!         print*,'tsurf_x,tsurf_w,tsurf,t1', ytsurf_th_x(j), ytsurf_th_w(j), ytsurf_th(j), yt(j,1)
2038         print*,'tsurf_x,t1x,tsurf_w,t1w,tsurf,t1,t1_ancien', &
2039 &               ytsurf_th_x(j), yt_x(j,1), ytsurf_th_w(j), yt_w(j,1), ytsurf_th(j), yt(j,1),t(j,1)
2040         print*,'qsatsurf,qsatsurf_x,qsatsurf_w', yqsatsurf(j), yqsatsurf_x(j), yqsatsurf_w(j)
2041         print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
2042        ENDDO
2043
2044        DO j=1,knon
2045         print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
2046 &             , y_flux_t1_x(j), y_flux_t1_w(j), y_flux_t1(j), y_flux_q1_x(j)*RLVTT, y_flux_q1_w(j)*RLVTT, yfluxlat(j), ywake_s(j)
2047         print*,'beta,ytsurf_new,yqsatsurf', ybeta(j), ytsurf_new(j), yqsatsurf(j)
2048         print*,'effusivity,facteur,cstar', effusivity, facteur,wake_cstar(j)
2049        ENDDO
2050       ENDIF  ! (prt_level >=10)
2051
2052!!! jyg le 07/02/2012
2053       ENDIF  ! (iflag_split .eq.1)
2054!!!
2055
2056!!! jyg le 07/02/2012
2057       IF (iflag_split .eq.0) THEN
2058!!!
2059!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2060        CALL climb_hq_up(knon, dtime, yt, yq, &
2061            y_flux_q1, y_flux_t1, ypaprs, ypplay, &
2062!!! jyg le 07/02/2012
2063            AcoefH, AcoefQ, BcoefH, BcoefQ, &
2064            CcoefH, CcoefQ, DcoefH, DcoefQ, &
2065            Kcoef_hq, gama_q, gama_h, &
2066!!!
2067            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))   
2068       ELSE  !(iflag_split .eq.0)
2069        CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
2070            y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
2071!!! nrlmd le 02/05/2011
2072            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
2073            CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
2074            Kcoef_hq_x, gama_q_x, gama_h_x, &
2075!!!
2076            y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:))   
2077!
2078       CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
2079            y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
2080!!! nrlmd le 02/05/2011
2081            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
2082            CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
2083            Kcoef_hq_w, gama_q_w, gama_h_w, &
2084!!!
2085            y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:))   
2086!!!
2087       ENDIF  ! (iflag_split .eq.0)
2088!!!
2089
2090!!! jyg le 07/02/2012
2091       IF (iflag_split .eq.0) THEN
2092!!!
2093!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2094        CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
2095!!! jyg le 07/02/2012
2096            AcoefU, AcoefV, BcoefU, BcoefV, &
2097            CcoefU, CcoefV, DcoefU, DcoefV, &
2098            Kcoef_m, &
2099!!!
2100            y_flux_u, y_flux_v, y_d_u, y_d_v)
2101     y_d_t_diss(:,:)=0.
2102     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2103        CALL yamada_c(knon,dtime,ypaprs,ypplay &
2104    &   ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
2105    &   ,iflag_pbl)
2106     ENDIF
2107!     print*,'yamada_c OK'
2108
2109       ELSE  !(iflag_split .eq.0)
2110        CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
2111!!! nrlmd le 02/05/2011
2112            AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
2113            CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
2114            Kcoef_m_x, &
2115!!!
2116            y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
2117!
2118     y_d_t_diss_x(:,:)=0.
2119     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2120        CALL yamada_c(knon,dtime,ypaprs,ypplay &
2121    &   ,yu_x,yv_x,yt_x,y_d_u_x,y_d_v_x,y_d_t_x,ycdragm_x,ytke_x,ycoefm_x,ycoefh_x &
2122        ,ycoefq_x,y_d_t_diss_x,yustar_x &
2123    &   ,iflag_pbl)
2124     ENDIF
2125!     print*,'yamada_c OK'
2126
2127        CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
2128!!! nrlmd le 02/05/2011
2129            AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
2130            CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
2131            Kcoef_m_w, &
2132!!!
2133            y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
2134!!!
2135     y_d_t_diss_w(:,:)=0.
2136     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
2137        CALL yamada_c(knon,dtime,ypaprs,ypplay &
2138    &   ,yu_w,yv_w,yt_w,y_d_u_w,y_d_v_w,y_d_t_w,ycdragm_w,ytke_w,ycoefm_w,ycoefh_w &
2139        ,ycoefq_w,y_d_t_diss_w,yustar_w &
2140    &   ,iflag_pbl)
2141     ENDIF
2142!     print*,'yamada_c OK'
2143!
2144        IF (prt_level >=10) THEN
2145         print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
2146               yfluxlat_x, yfluxlat_w
2147        ENDIF
2148!
2149       ENDIF  ! (iflag_split .eq.0)
2150!!!
2151
2152        DO j = 1, knon
2153          y_dflux_t(j) = y_dflux_t(j) * ypct(j)
2154          y_dflux_q(j) = y_dflux_q(j) * ypct(j)
2155        ENDDO
2156
2157!****************************************************************************************
2158! 13) Transform variables for output format :
2159!     - Decompress
2160!     - Multiply with pourcentage of current surface
2161!     - Cumulate in global variable
2162!
2163!****************************************************************************************
2164
2165
2166!!! jyg le 07/02/2012
2167       IF (iflag_split .eq.0) THEN
2168!!!
2169        DO k = 1, klev
2170           DO j = 1, knon
2171             i = ni(j)
2172             y_d_t_diss(j,k)  = y_d_t_diss(j,k) * ypct(j)
2173             y_d_t(j,k)  = y_d_t(j,k) * ypct(j)
2174             y_d_q(j,k)  = y_d_q(j,k) * ypct(j)
2175             y_d_u(j,k)  = y_d_u(j,k) * ypct(j)
2176             y_d_v(j,k)  = y_d_v(j,k) * ypct(j)
2177
2178             flux_t(i,k,nsrf) = y_flux_t(j,k)
2179             flux_q(i,k,nsrf) = y_flux_q(j,k)
2180             flux_u(i,k,nsrf) = y_flux_u(j,k)
2181             flux_v(i,k,nsrf) = y_flux_v(j,k)
2182
2183
2184           ENDDO
2185        ENDDO
2186
2187
2188       ELSE  !(iflag_split .eq.0)
2189
2190! Tendances hors poches
2191        DO k = 1, klev
2192          DO j = 1, knon
2193            i = ni(j)
2194            y_d_t_diss_x(j,k)  = y_d_t_diss_x(j,k) * ypct(j)
2195            y_d_t_x(j,k)  = y_d_t_x(j,k) * ypct(j)
2196            y_d_q_x(j,k)  = y_d_q_x(j,k) * ypct(j)
2197            y_d_u_x(j,k)  = y_d_u_x(j,k) * ypct(j)
2198            y_d_v_x(j,k)  = y_d_v_x(j,k) * ypct(j)
2199
2200            flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
2201            flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
2202            flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
2203            flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
2204          ENDDO
2205        ENDDO
2206
2207! Tendances dans les poches
2208        DO k = 1, klev
2209          DO j = 1, knon
2210            i = ni(j)
2211            y_d_t_diss_w(j,k)  = y_d_t_diss_w(j,k) * ypct(j)
2212            y_d_t_w(j,k)  = y_d_t_w(j,k) * ypct(j)
2213            y_d_q_w(j,k)  = y_d_q_w(j,k) * ypct(j)
2214            y_d_u_w(j,k)  = y_d_u_w(j,k) * ypct(j)
2215            y_d_v_w(j,k)  = y_d_v_w(j,k) * ypct(j)
2216
2217            flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
2218            flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
2219            flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
2220            flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
2221          ENDDO
2222        ENDDO
2223
2224! Flux, tendances et Tke moyenne dans la maille
2225        DO k = 1, klev
2226          DO j = 1, knon
2227            i = ni(j)
2228            flux_t(i,k,nsrf) = flux_t_x(i,k,nsrf)+ywake_s(j)*(flux_t_w(i,k,nsrf)-flux_t_x(i,k,nsrf))
2229            flux_q(i,k,nsrf) = flux_q_x(i,k,nsrf)+ywake_s(j)*(flux_q_w(i,k,nsrf)-flux_q_x(i,k,nsrf))
2230            flux_u(i,k,nsrf) = flux_u_x(i,k,nsrf)+ywake_s(j)*(flux_u_w(i,k,nsrf)-flux_u_x(i,k,nsrf))
2231            flux_v(i,k,nsrf) = flux_v_x(i,k,nsrf)+ywake_s(j)*(flux_v_w(i,k,nsrf)-flux_v_x(i,k,nsrf))
2232          ENDDO
2233        ENDDO
2234        DO j=1,knon
2235          yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
2236        ENDDO
2237        IF (prt_level >=10) THEN
2238          print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
2239                    nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
2240        ENDIF
2241
2242        DO k = 1, klev
2243          DO j = 1, knon
2244            y_d_t_diss(j,k) = y_d_t_diss_x(j,k)+ywake_s(j)*(y_d_t_diss_w(j,k) -y_d_t_diss_x(j,k))
2245            y_d_t(j,k) = y_d_t_x(j,k)+ywake_s(j)*(y_d_t_w(j,k) -y_d_t_x(j,k))
2246            y_d_q(j,k) = y_d_q_x(j,k)+ywake_s(j)*(y_d_q_w(j,k) -y_d_q_x(j,k))
2247            y_d_u(j,k) = y_d_u_x(j,k)+ywake_s(j)*(y_d_u_w(j,k) -y_d_u_x(j,k))
2248            y_d_v(j,k) = y_d_v_x(j,k)+ywake_s(j)*(y_d_v_w(j,k) -y_d_v_x(j,k))
2249          ENDDO
2250        ENDDO
2251
2252       ENDIF  ! (iflag_split .eq.0)
2253!!!
2254
2255!      print*,'Dans pbl OK1'
2256
2257!jyg<
2258!!       evap(:,nsrf) = - flux_q(:,1,nsrf)
2259!>jyg
2260       DO j = 1, knon
2261          i = ni(j)
2262          evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
2263          d_ts(i,nsrf) = y_d_ts(j)
2264!albedo SB >>>
2265          do k=1,nsw
2266          alb_dir(i,k,nsrf) = yalb_dir_new(j,k)
2267          alb_dif(i,k,nsrf) = yalb_dif_new(j,k)
2268          enddo
2269!albedo SB <<<
2270          snow(i,nsrf) = ysnow(j) 
2271          qsurf(i,nsrf) = yqsurf(j)
2272          z0m(i,nsrf) = yz0m(j)
2273          z0h(i,nsrf) = yz0h(j)
2274          fluxlat(i,nsrf) = yfluxlat(j)
2275          agesno(i,nsrf) = yagesno(j) 
2276          cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
2277          cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
2278          dflux_t(i) = dflux_t(i) + y_dflux_t(j)
2279          dflux_q(i) = dflux_q(i) + y_dflux_q(j)
2280       END DO
2281
2282!      print*,'Dans pbl OK2'
2283
2284!!! jyg le 07/02/2012
2285       IF (iflag_split .ge.1) THEN
2286!!!
2287!!! nrlmd le 02/05/2011
2288        DO j = 1, knon
2289          i = ni(j)
2290          fluxlat_x(i,nsrf) = yfluxlat_x(j)
2291          fluxlat_w(i,nsrf) = yfluxlat_w(j)
2292!!!
2293!!! nrlmd le 13/06/2011
2294          delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
2295          cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
2296          cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
2297          cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
2298          cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
2299          kh(i) = kh(i) + Kech_h(j)*ypct(j)
2300          kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
2301          kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
2302!!!
2303        END DO
2304!!!     
2305       ENDIF  ! (iflag_split .ge.1)
2306!!!
2307!!! nrlmd le 02/05/2011
2308!!jyg le 20/02/2011
2309!!        tke_x(:,:,nsrf)=0.
2310!!        tke_w(:,:,nsrf)=0.
2311!!jyg le 20/02/2011
2312!!        DO k = 1, klev+1
2313!!          DO j = 1, knon
2314!!            i = ni(j)
2315!!            wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2316!!            tke(i,k,nsrf)   = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2317!!          ENDDO
2318!!        ENDDO
2319!!jyg le 20/02/2011
2320!!        DO k = 1, klev+1
2321!!          DO j = 1, knon
2322!!            i = ni(j)
2323!!            tke(i,k,nsrf)=(1.-ywake_s(j))*tke_x(i,k,nsrf)+ywake_s(j)*tke_w(i,k,nsrf)
2324!!          ENDDO
2325!!        ENDDO
2326!!!
2327       IF (iflag_split .eq.0) THEN
2328        wake_dltke(:,:,nsrf) = 0.
2329        DO k = 1, klev
2330           DO j = 1, knon
2331              i = ni(j)
2332!jyg<
2333!!              tke(i,k,nsrf)    = ytke(j,k)
2334!!              tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
2335              tke_x(i,k,nsrf)    = ytke(j,k)
2336              tke_x(i,k,is_ave) = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
2337!>jyg
2338           END DO
2339        END DO
2340
2341       ELSE  ! (iflag_split .eq.0)
2342        DO k = 1, klev
2343          DO j = 1, knon
2344            i = ni(j)
2345            wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
2346!jyg<
2347!!            tke(i,k,nsrf)   = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
2348!!            tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
2349            tke_x(i,k,nsrf)   = ytke_x(j,k)
2350            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)
2351            wake_dltke(i,k,is_ave)   = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
2352
2353!>jyg
2354          ENDDO
2355        ENDDO
2356       ENDIF  ! (iflag_split .eq.0)
2357!!!
2358       DO k = 2, klev
2359          DO j = 1, knon
2360             i = ni(j)
2361             zcoefh(i,k,nsrf) = ycoefh(j,k)
2362             zcoefm(i,k,nsrf) = ycoefm(j,k)
2363             zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
2364             zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
2365          END DO
2366       END DO
2367
2368!      print*,'Dans pbl OK3'
2369
2370       IF ( nsrf .EQ. is_ter ) THEN
2371          DO j = 1, knon
2372             i = ni(j)
2373             qsol(i) = yqsol(j)
2374          END DO
2375       END IF
2376       
2377!jyg<
2378!!       ftsoil(:,:,nsrf) = 0.
2379!>jyg
2380       DO k = 1, nsoilmx
2381          DO j = 1, knon
2382             i = ni(j)
2383             ftsoil(i, k, nsrf) = ytsoil(j,k)
2384          END DO
2385       END DO
2386       
2387!!! jyg le 07/02/2012
2388       IF (iflag_split .ge.1) THEN
2389!!!
2390!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
2391        DO k = 1, klev
2392          DO j = 1, knon
2393           i = ni(j)
2394           d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
2395           d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
2396           d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
2397           d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
2398           d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
2399!
2400           d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
2401           d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
2402           d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
2403           d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
2404           d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
2405!
2406!!           d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
2407!!           d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
2408          END DO
2409        END DO
2410!!!
2411       ENDIF  ! (iflag_split .ge.1)
2412!!!
2413       
2414       DO k = 1, klev
2415          DO j = 1, knon
2416             i = ni(j)
2417             d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
2418             d_t(i,k) = d_t(i,k) + y_d_t(j,k)
2419             d_q(i,k) = d_q(i,k) + y_d_q(j,k)
2420             d_u(i,k) = d_u(i,k) + y_d_u(j,k)
2421             d_v(i,k) = d_v(i,k) + y_d_v(j,k)
2422          END DO
2423       END DO
2424
2425!      print*,'Dans pbl OK4'
2426
2427       IF (prt_level >=10) THEN
2428         print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
2429          d_t_w(:,1), d_t_x(:,1), d_t(:,1)
2430       ENDIF
2431
2432!****************************************************************************************
2433! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
2434!     Call HBTM
2435!
2436!****************************************************************************************
2437!!!
2438!
2439#undef T2m     
2440#define T2m     
2441#ifdef T2m
2442! Calculations of diagnostic t,q at 2m and u, v at 10m
2443
2444!      print*,'Dans pbl OK41'
2445!      print*,'tair1,yt(:,1),y_d_t(:,1)'
2446!      print*, tair1,yt(:,1),y_d_t(:,1)
2447!!! jyg le 07/02/2012
2448       IF (iflag_split .eq.0) THEN
2449        DO j=1, knon
2450          uzon(j) = yu(j,1) + y_d_u(j,1)
2451          vmer(j) = yv(j,1) + y_d_v(j,1)
2452          tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
2453          qair1(j) = yq(j,1) + y_d_q(j,1)
2454          zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2455               * (ypaprs(j,1)-ypplay(j,1))
2456          tairsol(j) = yts(j) + y_d_ts(j)
2457          qairsol(j) = yqsurf(j)
2458        END DO
2459       ELSE  ! (iflag_split .eq.0)
2460        DO j=1, knon
2461          uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
2462          vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
2463          tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
2464          qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
2465          zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2466               * (ypaprs(j,1)-ypplay(j,1))
2467          tairsol(j) = yts(j) + y_d_ts(j)
2468          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
2469          qairsol(j) = yqsurf(j)
2470        END DO
2471        DO j=1, knon
2472          uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
2473          vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
2474          tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
2475          qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
2476          zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
2477               * (ypaprs(j,1)-ypplay(j,1))
2478          tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
2479          qairsol(j) = yqsurf(j)
2480        END DO
2481!!!     
2482       ENDIF  ! (iflag_split .eq.0)
2483!!!
2484       DO j=1, knon
2485          i = ni(j)
2486          rugo1(j) = yz0m(j)
2487          IF(nsrf.EQ.is_oce) THEN
2488             rugo1(j) = z0m(i,nsrf)
2489          ENDIF
2490          psfce(j)=ypaprs(j,1)
2491          patm(j)=ypplay(j,1)
2492       END DO
2493       
2494!      print*,'Dans pbl OK42A'
2495!      print*,'tair1,yt(:,1),y_d_t(:,1)'
2496!      print*, tair1,yt(:,1),y_d_t(:,1)
2497
2498! Calculate the temperature et relative humidity at 2m and the wind at 10m
2499!!! jyg le 07/02/2012
2500       IF (iflag_split .eq.0) THEN
2501        CALL stdlevvar(klon, knon, nsrf, zxli, &
2502            uzon, vmer, tair1, qair1, zgeo1, &
2503            tairsol, qairsol, rugo1, rugo1, psfce, patm, &
2504            yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
2505       ELSE  !(iflag_split .eq.0)
2506        CALL stdlevvar(klon, knon, nsrf, zxli, &
2507            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
2508            tairsol_x, qairsol, rugo1, rugo1, psfce, patm, &
2509            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x)
2510        CALL stdlevvar(klon, knon, nsrf, zxli, &
2511            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
2512            tairsol_w, qairsol, rugo1, rugo1, psfce, patm, &
2513            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w)
2514!!!
2515       ENDIF  ! (iflag_split .eq.0)
2516!!!
2517!!! jyg le 07/02/2012
2518       IF (iflag_split .eq.0) THEN
2519        DO j=1, knon
2520          i = ni(j)
2521          t2m(i,nsrf)=yt2m(j)
2522          q2m(i,nsrf)=yq2m(j)
2523     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2524          ustar(i,nsrf)=yustar(j)
2525          u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
2526          v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
2527        END DO
2528       ELSE  !(iflag_split .eq.0)
2529        DO j=1, knon
2530          i = ni(j)
2531          t2m_x(i,nsrf)=yt2m_x(j)
2532          q2m_x(i,nsrf)=yq2m_x(j)
2533     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2534          ustar_x(i,nsrf)=yustar_x(j)
2535          u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
2536          v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
2537        END DO
2538        DO j=1, knon
2539          i = ni(j)
2540          t2m_w(i,nsrf)=yt2m_w(j)
2541          q2m_w(i,nsrf)=yq2m_w(j)
2542     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
2543          ustar_w(i,nsrf)=yustar_w(j)
2544          u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)**2+vmer_w(j)**2)
2545          v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)**2+vmer_w(j)**2)
2546!
2547          ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
2548          u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
2549          v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
2550        END DO
2551!!!
2552       ENDIF  ! (iflag_split .eq.0)
2553!!!
2554
2555!      print*,'Dans pbl OK43'
2556!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
2557!IM Ajoute dependance type surface
2558       IF (thermcep) THEN
2559!!! jyg le 07/02/2012
2560       IF (iflag_split .eq.0) THEN
2561          DO j = 1, knon
2562             i=ni(j)
2563             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
2564             zx_qs1  = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
2565             zx_qs1  = MIN(0.5,zx_qs1)
2566             zcor1   = 1./(1.-RETV*zx_qs1)
2567             zx_qs1  = zx_qs1*zcor1
2568             
2569             rh2m(i)   = rh2m(i)   + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
2570             qsat2m(i) = qsat2m(i) + zx_qs1  * pctsrf(i,nsrf)
2571          END DO
2572       ELSE  ! (iflag_split .eq.0)
2573          DO j = 1, knon
2574             i=ni(j)
2575             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
2576             zx_qs1  = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
2577             zx_qs1  = MIN(0.5,zx_qs1)
2578             zcor1   = 1./(1.-RETV*zx_qs1)
2579             zx_qs1  = zx_qs1*zcor1
2580             
2581             rh2m_x(i)   = rh2m_x(i)   + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
2582             qsat2m_x(i) = qsat2m_x(i) + zx_qs1  * pctsrf(i,nsrf)
2583          END DO
2584          DO j = 1, knon
2585             i=ni(j)
2586             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
2587             zx_qs1  = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
2588             zx_qs1  = MIN(0.5,zx_qs1)
2589             zcor1   = 1./(1.-RETV*zx_qs1)
2590             zx_qs1  = zx_qs1*zcor1
2591             
2592             rh2m_w(i)   = rh2m_w(i)   + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
2593             qsat2m_w(i) = qsat2m_w(i) + zx_qs1  * pctsrf(i,nsrf)
2594          END DO
2595!!!     
2596       ENDIF  ! (iflag_split .eq.0)
2597!!!
2598       END IF
2599!
2600       IF (prt_level >=10) THEN
2601         print *, 'T2m, q2m, RH2m ', &
2602          t2m, q2m, rh2m
2603       ENDIF
2604
2605!   print*,'OK pbl 5'
2606!
2607!!! jyg le 07/02/2012
2608       IF (iflag_split .eq.0) THEN
2609        CALL hbtm(knon, ypaprs, ypplay, &
2610            yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
2611            y_flux_t,y_flux_q,yu,yv,yt,yq, &
2612            ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
2613            ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
2614          IF (prt_level >=10) THEN
2615       print *,' Arg. de HBTM: yt2m ',yt2m
2616       print *,' Arg. de HBTM: yt10m ',yt10m
2617       print *,' Arg. de HBTM: yq2m ',yq2m
2618       print *,' Arg. de HBTM: yq10m ',yq10m
2619       print *,' Arg. de HBTM: yustar ',yustar
2620       print *,' Arg. de HBTM: y_flux_t ',y_flux_t
2621       print *,' Arg. de HBTM: y_flux_q ',y_flux_q
2622       print *,' Arg. de HBTM: yu ',yu
2623       print *,' Arg. de HBTM: yv ',yv
2624       print *,' Arg. de HBTM: yt ',yt
2625       print *,' Arg. de HBTM: yq ',yq
2626          ENDIF
2627       ELSE  ! (iflag_split .eq.0)
2628        CALL HBTM(knon, ypaprs, ypplay, &
2629            yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
2630            y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
2631            ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
2632            ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
2633          IF (prt_level >=10) THEN
2634       print *,' Arg. de HBTM: yt2m_x ',yt2m_x
2635       print *,' Arg. de HBTM: yt10m_x ',yt10m_x
2636       print *,' Arg. de HBTM: yq2m_x ',yq2m_x
2637       print *,' Arg. de HBTM: yq10m_x ',yq10m_x
2638       print *,' Arg. de HBTM: yustar_x ',yustar_x
2639       print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x
2640       print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x
2641       print *,' Arg. de HBTM: yu_x ',yu_x
2642       print *,' Arg. de HBTM: yv_x ',yv_x
2643       print *,' Arg. de HBTM: yt_x ',yt_x
2644       print *,' Arg. de HBTM: yq_x ',yq_x
2645          ENDIF
2646        CALL HBTM(knon, ypaprs, ypplay, &
2647            yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
2648            y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
2649            ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
2650            ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
2651!!!     
2652       ENDIF  ! (iflag_split .eq.0)
2653!!!
2654       
2655!!! jyg le 07/02/2012
2656       IF (iflag_split .eq.0) THEN
2657!!!
2658        DO j=1, knon
2659          i = ni(j)
2660          pblh(i,nsrf)   = ypblh(j)
2661          wstar(i,nsrf)  = ywstar(j)
2662          plcl(i,nsrf)   = ylcl(j)
2663          capCL(i,nsrf)  = ycapCL(j)
2664          oliqCL(i,nsrf) = yoliqCL(j)
2665          cteiCL(i,nsrf) = ycteiCL(j)
2666          pblT(i,nsrf)   = ypblT(j)
2667          therm(i,nsrf)  = ytherm(j)
2668          trmb1(i,nsrf)  = ytrmb1(j)
2669          trmb2(i,nsrf)  = ytrmb2(j)
2670          trmb3(i,nsrf)  = ytrmb3(j)
2671        END DO
2672        IF (prt_level >=10) THEN
2673          print *, 'After HBTM: pblh ', pblh
2674          print *, 'After HBTM: plcl ', plcl
2675          print *, 'After HBTM: cteiCL ', cteiCL
2676        ENDIF
2677       ELSE  !(iflag_split .eq.0)
2678        DO j=1, knon
2679          i = ni(j)
2680          pblh_x(i,nsrf)   = ypblh_x(j)
2681          wstar_x(i,nsrf)  = ywstar_x(j)
2682          plcl_x(i,nsrf)   = ylcl_x(j)
2683          capCL_x(i,nsrf)  = ycapCL_x(j)
2684          oliqCL_x(i,nsrf) = yoliqCL_x(j)
2685          cteiCL_x(i,nsrf) = ycteiCL_x(j)
2686          pblT_x(i,nsrf)   = ypblT_x(j)
2687          therm_x(i,nsrf)  = ytherm_x(j)
2688          trmb1_x(i,nsrf)  = ytrmb1_x(j)
2689          trmb2_x(i,nsrf)  = ytrmb2_x(j)
2690          trmb3_x(i,nsrf)  = ytrmb3_x(j)
2691        END DO
2692        IF (prt_level >=10) THEN
2693          print *, 'After HBTM: pblh_x ', pblh_x
2694          print *, 'After HBTM: plcl_x ', plcl_x
2695          print *, 'After HBTM: cteiCL_x ', cteiCL_x
2696        ENDIF
2697        DO j=1, knon
2698          i = ni(j)
2699          pblh_w(i,nsrf)   = ypblh_w(j)
2700          wstar_w(i,nsrf)  = ywstar_w(j)
2701          plcl_w(i,nsrf)   = ylcl_w(j)
2702          capCL_w(i,nsrf)  = ycapCL_w(j)
2703          oliqCL_w(i,nsrf) = yoliqCL_w(j)
2704          cteiCL_w(i,nsrf) = ycteiCL_w(j)
2705          pblT_w(i,nsrf)   = ypblT_w(j)
2706          therm_w(i,nsrf)  = ytherm_w(j)
2707          trmb1_w(i,nsrf)  = ytrmb1_w(j)
2708          trmb2_w(i,nsrf)  = ytrmb2_w(j)
2709          trmb3_w(i,nsrf)  = ytrmb3_w(j)
2710        END DO
2711        IF (prt_level >=10) THEN
2712          print *, 'After HBTM: pblh_w ', pblh_w
2713          print *, 'After HBTM: plcl_w ', plcl_w
2714          print *, 'After HBTM: cteiCL_w ', cteiCL_w
2715        ENDIF
2716!!!
2717       ENDIF  ! (iflag_split .eq.0)
2718!!!
2719
2720!   print*,'OK pbl 6'
2721#else
2722! T2m not defined
2723! No calculation
2724       PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
2725#endif
2726
2727!****************************************************************************************
2728! 15) End of loop over different surfaces
2729!
2730!****************************************************************************************
2731    END DO loop_nbsrf
2732
2733!****************************************************************************************
2734! 16) Calculate the mean value over all sub-surfaces for some variables
2735!
2736!****************************************************************************************
2737   
2738    z0m(:,nbsrf+1) = 0.0
2739    z0h(:,nbsrf+1) = 0.0
2740    DO nsrf = 1, nbsrf
2741       DO i = 1, klon
2742          z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
2743          z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
2744       ENDDO
2745    ENDDO
2746
2747!   print*,'OK pbl 7'
2748    zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
2749    zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
2750    zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
2751    zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
2752    zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
2753    zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
2754
2755!!! jyg le 07/02/2012
2756       IF (iflag_split .ge.1) THEN
2757!!!
2758!!! nrlmd & jyg les 02/05/2011, 05/02/2012
2759
2760        DO nsrf = 1, nbsrf
2761          DO k = 1, klev
2762            DO i = 1, klon
2763              zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
2764              zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
2765              zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
2766              zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
2767!
2768              zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
2769              zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
2770              zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
2771              zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
2772            END DO
2773          END DO
2774        END DO
2775
2776    DO i = 1, klon
2777      zxsens_x(i) = - zxfluxt_x(i,1)
2778      zxsens_w(i) = - zxfluxt_w(i,1)
2779    END DO
2780!!!
2781       ENDIF  ! (iflag_split .ge.1)
2782!!!
2783
2784    DO nsrf = 1, nbsrf
2785       DO k = 1, klev
2786          DO i = 1, klon
2787             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
2788             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
2789             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
2790             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
2791          END DO
2792       END DO
2793    END DO
2794
2795    DO i = 1, klon
2796       zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
2797       zxevap(i)     = - zxfluxq(i,1) ! flux d'evaporation au sol
2798       fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
2799    ENDDO
2800!!!
2801   
2802!
2803! Incrementer la temperature du sol
2804!
2805    zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
2806    zt2m(:) = 0.0    ; zq2m(:) = 0.0
2807    zustar(:)=0.0 ; zu10m(:) = 0.0   ; zv10m(:) = 0.0
2808    s_pblh(:) = 0.0  ; s_plcl(:) = 0.0
2809!!! jyg le 07/02/2012
2810     s_pblh_x(:) = 0.0  ; s_plcl_x(:) = 0.0
2811     s_pblh_w(:) = 0.0  ; s_plcl_w(:) = 0.0
2812!!!
2813    s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
2814    s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
2815    s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
2816    s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
2817    wstar(:,is_ave)=0.
2818   
2819!   print*,'OK pbl 9'
2820   
2821!!! nrlmd le 02/05/2011
2822    zxfluxlat_x(:) = 0.0  ;  zxfluxlat_w(:) = 0.0
2823!!!
2824   
2825    DO nsrf = 1, nbsrf
2826       DO i = 1, klon         
2827          ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
2828         
2829          wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
2830               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
2831          wfbilo(i,nsrf) = (evap(i,nsrf)-(rain_f(i)+snow_f(i)))*pctsrf(i,nsrf)
2832          wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
2833          wfrain(i,nsrf) = rain_f(i)*pctsrf(i,nsrf)
2834          wfsnow(i,nsrf) = snow_f(i)*pctsrf(i,nsrf)
2835
2836          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf(i,nsrf)
2837          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
2838       END DO
2839    END DO
2840!
2841!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
2842   IF (iflag_order2_sollw == 1) THEN
2843    meansqT(:) = 0. ! as working buffer
2844    DO nsrf = 1, nbsrf
2845     DO i = 1, klon
2846      meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))**2 *pctsrf(i,nsrf)
2847     END DO
2848    END DO
2849    zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
2850   ENDIF   ! iflag_order2_sollw == 1
2851!>al1
2852         
2853!!! jyg le 07/02/2012
2854       IF (iflag_split .eq.0) THEN
2855        DO nsrf = 1, nbsrf
2856         DO i = 1, klon         
2857          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
2858          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
2859          zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
2860          wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
2861          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
2862          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
2863
2864          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf(i,nsrf)
2865          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf(i,nsrf)
2866          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf(i,nsrf)
2867          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
2868          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
2869          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf(i,nsrf)
2870          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf(i,nsrf)
2871          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf(i,nsrf)
2872          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf(i,nsrf)
2873          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf(i,nsrf)
2874         END DO
2875        END DO
2876       ELSE  !(iflag_split .eq.0)
2877        DO nsrf = 1, nbsrf
2878         DO i = 1, klon         
2879!!! nrlmd le 02/05/2011
2880          zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
2881          zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
2882!!!
2883!!! jyg le 08/02/2012
2884!!  Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
2885!!  pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
2886!!  pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
2887!!  pour les autres variables, on sort les valeurs de la region (x).
2888          zt2m(i)  = zt2m(i)  + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf)
2889          zq2m(i)  = zq2m(i)  + q2m_x(i,nsrf)  * pctsrf(i,nsrf)
2890          zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
2891          wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
2892          zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
2893          zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
2894!
2895          s_pblh(i)     = s_pblh(i)     + pblh_x(i,nsrf)  * pctsrf(i,nsrf)
2896          s_pblh_x(i)   = s_pblh_x(i)   + pblh_x(i,nsrf)  * pctsrf(i,nsrf)
2897          s_pblh_w(i)   = s_pblh_w(i)   + pblh_w(i,nsrf)  * pctsrf(i,nsrf)
2898!
2899          s_plcl(i)     = s_plcl(i)     + plcl_x(i,nsrf)  * pctsrf(i,nsrf)
2900          s_plcl_x(i)   = s_plcl_x(i)   + plcl_x(i,nsrf)  * pctsrf(i,nsrf)
2901          s_plcl_w(i)   = s_plcl_w(i)   + plcl_w(i,nsrf)  * pctsrf(i,nsrf)
2902!
2903          s_capCL(i)  = s_capCL(i)  + capCL_x(i,nsrf) * pctsrf(i,nsrf)
2904          s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
2905          s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
2906          s_pblT(i)   = s_pblT(i)   + pblT_x(i,nsrf)  * pctsrf(i,nsrf)
2907          s_therm(i)  = s_therm(i)  + therm_x(i,nsrf) * pctsrf(i,nsrf)
2908          s_trmb1(i)  = s_trmb1(i)  + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
2909          s_trmb2(i)  = s_trmb2(i)  + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
2910          s_trmb3(i)  = s_trmb3(i)  + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
2911         END DO
2912        END DO
2913        DO i = 1, klon         
2914          qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
2915        END DO
2916!!!
2917       ENDIF  ! (iflag_split .eq.0)
2918!!!
2919
2920    IF (check) THEN
2921       amn=MIN(ts(1,is_ter),1000.)
2922       amx=MAX(ts(1,is_ter),-1000.)
2923       DO i=2, klon
2924          amn=MIN(ts(i,is_ter),amn)
2925          amx=MAX(ts(i,is_ter),amx)
2926       ENDDO
2927       PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
2928    ENDIF
2929
2930!jg ?
2931!!$!
2932!!$! If a sub-surface does not exsist for a grid point, the mean value for all
2933!!$! sub-surfaces is distributed.
2934!!$!
2935!!$    DO nsrf = 1, nbsrf
2936!!$       DO i = 1, klon
2937!!$          IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
2938!!$             ts(i,nsrf)     = zxtsol(i)
2939!!$             t2m(i,nsrf)    = zt2m(i)
2940!!$             q2m(i,nsrf)    = zq2m(i)
2941!!$             u10m(i,nsrf)   = zu10m(i)
2942!!$             v10m(i,nsrf)   = zv10m(i)
2943!!$
2944!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
2945!!$             pblh(i,nsrf)   = s_pblh(i)
2946!!$             plcl(i,nsrf)   = s_plcl(i)
2947!!$             capCL(i,nsrf)  = s_capCL(i)
2948!!$             oliqCL(i,nsrf) = s_oliqCL(i)
2949!!$             cteiCL(i,nsrf) = s_cteiCL(i)
2950!!$             pblT(i,nsrf)   = s_pblT(i)
2951!!$             therm(i,nsrf)  = s_therm(i)
2952!!$             trmb1(i,nsrf)  = s_trmb1(i)
2953!!$             trmb2(i,nsrf)  = s_trmb2(i)
2954!!$             trmb3(i,nsrf)  = s_trmb3(i)
2955!!$          ENDIF
2956!!$       ENDDO
2957!!$    ENDDO
2958
2959
2960    DO i = 1, klon
2961       fder(i) = - 4.0*RSIGMA*zxtsol(i)**3
2962    ENDDO
2963   
2964    zxqsurf(:) = 0.0
2965    zxsnow(:)  = 0.0
2966    DO nsrf = 1, nbsrf
2967       DO i = 1, klon
2968          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
2969          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
2970       END DO
2971    END DO
2972
2973! Premier niveau de vent sortie dans physiq.F
2974    zu1(:) = u(:,1)
2975    zv1(:) = v(:,1)
2976
2977
2978  END SUBROUTINE pbl_surface
2979!
2980!****************************************************************************************
2981!
2982  SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
2983
2984    USE indice_sol_mod
2985
2986    INCLUDE "dimsoil.h"
2987
2988! Ouput variables
2989!****************************************************************************************
2990    REAL, DIMENSION(klon), INTENT(OUT)                 :: fder_rst
2991    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: snow_rst
2992    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
2993    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
2994
2995 
2996!****************************************************************************************
2997! Return module variables for writing to restart file
2998!
2999!****************************************************************************************   
3000    fder_rst(:)       = fder(:)
3001    snow_rst(:,:)     = snow(:,:)
3002    qsurf_rst(:,:)    = qsurf(:,:)
3003    ftsoil_rst(:,:,:) = ftsoil(:,:,:)
3004
3005!****************************************************************************************
3006! Deallocate module variables
3007!
3008!****************************************************************************************
3009!   DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
3010    IF (ALLOCATED(fder)) DEALLOCATE(fder)
3011    IF (ALLOCATED(snow)) DEALLOCATE(snow)
3012    IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
3013    IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
3014
3015  END SUBROUTINE pbl_surface_final
3016
3017!****************************************************************************************
3018!
3019
3020!albedo SB >>>
3021SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
3022     evap, z0m, z0h, agesno,                                  &
3023     tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke) 
3024!albedo SB <<<
3025    ! Give default values where new fraction has appread
3026
3027    USE indice_sol_mod
3028
3029    INCLUDE "dimsoil.h"
3030    INCLUDE "clesphys.h"
3031    INCLUDE "compbl.h"
3032
3033! Input variables
3034!****************************************************************************************
3035    INTEGER, INTENT(IN)                     :: itime
3036    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
3037
3038! InOutput variables
3039!****************************************************************************************
3040    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: tsurf
3041!albedo SB >>>
3042    REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT)       :: alb_dir, alb_dif
3043    INTEGER :: k
3044!albedo SB <<<
3045    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: ustar,u10m, v10m
3046    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: evap, agesno
3047    REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT)        :: z0m,z0h
3048    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
3049
3050! Local variables
3051!****************************************************************************************
3052    INTEGER           :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
3053    CHARACTER(len=80) :: abort_message
3054    CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
3055    INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
3056!
3057! All at once !!
3058!****************************************************************************************
3059   
3060    DO nsrf = 1, nbsrf
3061       ! First decide complement sub-surfaces
3062       SELECT CASE (nsrf)
3063       CASE(is_oce)
3064          nsrf_comp1=is_sic
3065          nsrf_comp2=is_ter
3066          nsrf_comp3=is_lic
3067       CASE(is_sic)
3068          nsrf_comp1=is_oce
3069          nsrf_comp2=is_ter
3070          nsrf_comp3=is_lic
3071       CASE(is_ter)
3072          nsrf_comp1=is_lic
3073          nsrf_comp2=is_oce
3074          nsrf_comp3=is_sic
3075       CASE(is_lic)
3076          nsrf_comp1=is_ter
3077          nsrf_comp2=is_oce
3078          nsrf_comp3=is_sic
3079       END SELECT
3080
3081       ! Initialize all new fractions
3082       DO i=1, klon
3083          IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
3084             
3085             IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
3086                ! Use the complement sub-surface, keeping the continents unchanged
3087                qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
3088                evap(i,nsrf)  = evap(i,nsrf_comp1)
3089                z0m(i,nsrf) = z0m(i,nsrf_comp1)
3090                z0h(i,nsrf) = z0h(i,nsrf_comp1)
3091                tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
3092!albedo SB >>>
3093                DO k=1,nsw
3094                 alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
3095                 alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
3096                ENDDO
3097!albedo SB <<<
3098                ustar(i,nsrf)  = ustar(i,nsrf_comp1)
3099                u10m(i,nsrf)  = u10m(i,nsrf_comp1)
3100                v10m(i,nsrf)  = v10m(i,nsrf_comp1)
3101                if (iflag_pbl > 1) then
3102                 tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
3103                endif
3104                mfois(nsrf) = mfois(nsrf) + 1
3105                ! F. Codron sensible default values for ocean and sea ice
3106                IF (nsrf.EQ.is_oce) THEN
3107                    tsurf(i,nsrf) = 271.35 ! Freezing sea water
3108                    DO k=1,nsw
3109                      alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
3110                      alb_dif(i,k,nsrf) = 0.06
3111                    END DO
3112                ELSE IF (nsrf.EQ.is_sic) THEN
3113                    tsurf(i,nsrf) = 273.15 ! Melting ice
3114                    DO k=1,nsw
3115                      alb_dir(i,k,nsrf) = 0.3 ! thin ice
3116                      alb_dif(i,k,nsrf) = 0.3
3117                    END DO
3118                END IF
3119                ! F. Codron
3120             ELSE
3121                ! The continents have changed. The new fraction receives the mean sum of the existent fractions
3122                qsurf(i,nsrf) = qsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3123                evap(i,nsrf)  = evap(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
3124                z0m(i,nsrf) = z0m(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3125                z0h(i,nsrf) = z0h(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3126                tsurf(i,nsrf) = tsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3127!albedo SB >>>
3128                DO k=1,nsw
3129                 alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3130                                        alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3131                 alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
3132                                        alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3133                ENDDO
3134!albedo SB <<<
3135                ustar(i,nsrf)  = ustar(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
3136                u10m(i,nsrf)  = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
3137                v10m(i,nsrf)  = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
3138                if (iflag_pbl > 1) then
3139                 tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
3140                endif
3141           
3142                ! Security abort. This option has never been tested. To test, comment the following line.
3143!                abort_message='The fraction of the continents have changed!'
3144!                CALL abort_physic(modname,abort_message,1)
3145                nfois(nsrf) = nfois(nsrf) + 1
3146             END IF
3147             snow(i,nsrf)     = 0.
3148             agesno(i,nsrf)   = 0.
3149             ftsoil(i,:,nsrf) = tsurf(i,nsrf)
3150          ELSE
3151             pfois(nsrf) = pfois(nsrf)+ 1
3152          END IF
3153       END DO
3154       
3155    END DO
3156
3157  END SUBROUTINE pbl_surface_newfrac
3158
3159
3160!****************************************************************************************
3161
3162
3163END MODULE pbl_surface_mod
3164
Note: See TracBrowser for help on using the repository browser.