source: LMDZ6/branches/LMDZ-INCA-Dyn/libf/phylmd/pbl_surface_mod.F90 @ 5456

Last change on this file since 5456 was 3784, checked in by evignon, 4 years ago

correction v3783 pour convergence et compilation ancienne physique, Etienne aide par Ehouarn

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