source: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 @ 3761

Last change on this file since 3761 was 3756, checked in by oboucher, 4 years ago

The fraction of SW down radiation that is diffuse is extracted from radlwsw.F90 and passed on to pbl_surface for inclusion in fields_out for transfer to ORCHIDEE if it requested in input file coupling_fields.def. Hence there is nothing automatic here. The fraction of SW down radiation is also added as a diagnostic. To satisfy the case when this new quantity is used in ORCHIDEE, then it is added in the restart (and read in case it is there).

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