source: LMDZ6/branches/contrails/libf/phylmd/pbl_surface_mod.F90 @ 5450

Last change on this file since 5450 was 5310, checked in by abarral, 7 weeks ago

unify abort_gcm
rename wxios -> wxios_mod

  • 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: 174.3 KB
Line 
1!
2! $Id: pbl_surface_mod.F90 5310 2024-11-01 12:05:47Z 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, landice_opt
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_qbs_mod,       ONLY : climb_qbs_down, climb_qbs_up
24  USE climb_wind_mod,      ONLY : climb_wind_down, climb_wind_up
25  USE coef_diff_turb_mod,  ONLY : coef_diff_turb
26  USE lmdz_call_atke,      ONLY : call_atke
27  USE ioipsl_getin_p_mod,  ONLY : getin_p
28  USE cdrag_mod
29  USE stdlevvar_mod
30  USE wx_pbl_var_mod,      ONLY : wx_pbl_init, wx_pbl_final, &
31                                  wx_pbl_prelim_0, wx_pbl_prelim_beta
32  USE wx_pbl_mod,          ONLY : wx_pbl0_merge, wx_pbl_split, wx_pbl_dts_merge, &
33                                  wx_pbl_check, wx_pbl_dts_check, wx_evappot
34  use config_ocean_skin_m, only: activate_ocean_skin
35#ifdef ISO
36  USE infotrac_phy, ONLY: niso,ntraciso=>ntiso   
37#endif
38
39  IMPLICIT NONE
40
41! Declaration of variables saved in restart file
42  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
43  !$OMP THREADPRIVATE(fder)
44  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE    :: snow   ! snow at surface
45  !$OMP THREADPRIVATE(snow)
46  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface
47  !$OMP THREADPRIVATE(qsurf)
48  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE          :: ftsoil ! soil temperature
49  !$OMP THREADPRIVATE(ftsoil)
50  REAL, ALLOCATABLE, DIMENSION(:), SAVE              :: ydTs0, ydqs0 
51                                                     ! nul forced temperature and humidity differences
52  !$OMP THREADPRIVATE(ydTs0, ydqs0)
53
54#ifdef ISO
55  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE   :: xtsnow   ! snow at surface
56  !$OMP THREADPRIVATE(xtsnow)
57  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: Rland_ice   ! snow at surface
58  !$OMP THREADPRIVATE(Rland_ice) 
59  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: Roce   ! snow at surface
60  !$OMP THREADPRIVATE(Roce) 
61#endif
62
63  INTEGER, SAVE :: iflag_pbl_surface_t2m_bug
64  !$OMP THREADPRIVATE(iflag_pbl_surface_t2m_bug)
65  INTEGER, SAVE :: iflag_new_t2mq2m
66  !$OMP THREADPRIVATE(iflag_new_t2mq2m)
67  LOGICAL, SAVE :: ok_bug_zg_wk_pbl
68  !$OMP THREADPRIVATE(ok_bug_zg_wk_pbl)
69
70!FC
71!  integer, save :: iflag_frein
72!  !$OMP THREADPRIVATE(iflag_frein)
73
74CONTAINS
75!
76!****************************************************************************************
77!
78  SUBROUTINE pbl_surface_init(fder_rst, snow_rst, qsurf_rst, ftsoil_rst)
79
80! This routine should be called after the restart file has been read.
81! This routine initialize the restart variables and does some validation tests
82! for the index of the different surfaces and tests the choice of type of ocean.
83
84    USE indice_sol_mod
85    USE print_control_mod, ONLY: lunout
86    USE ioipsl_getin_p_mod, ONLY : getin_p
87    USE dimsoil_mod_h, ONLY: nsoilmx
88    IMPLICIT NONE
89 
90! Input variables
91!****************************************************************************************
92    REAL, DIMENSION(klon), INTENT(IN)                 :: fder_rst
93    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: snow_rst
94    REAL, DIMENSION(klon, nbsrf), INTENT(IN)          :: qsurf_rst
95    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
96 
97! Local variables
98!****************************************************************************************
99    INTEGER                       :: ierr
100    CHARACTER(len=80)             :: abort_message
101    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
102   
103!****************************************************************************************
104! Allocate and initialize module variables with fields read from restart file.
105!
106!****************************************************************************************   
107    ALLOCATE(fder(klon), stat=ierr)
108    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
109
110    ALLOCATE(snow(klon,nbsrf), stat=ierr)
111    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
112
113    ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
114    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
115
116    ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
117    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
118
119    ALLOCATE(ydTs0(klon), stat=ierr)
120    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
121
122    ALLOCATE(ydqs0(klon), stat=ierr)
123    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
124
125    fder(:)       = fder_rst(:)
126    snow(:,:)     = snow_rst(:,:)
127    qsurf(:,:)    = qsurf_rst(:,:)
128    ftsoil(:,:,:) = ftsoil_rst(:,:,:)
129    ydTs0(:) = 0.
130    ydqs0(:) = 0.
131
132!****************************************************************************************
133! Test for sub-surface indices
134!
135!****************************************************************************************
136    IF (is_ter /= 1) THEN
137      WRITE(lunout,*)" *** Warning ***"
138      WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
139      WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
140      abort_message="voir ci-dessus"
141      CALL abort_physic(modname,abort_message,1)
142    ENDIF
143
144    IF ( is_oce > is_sic ) THEN
145      WRITE(lunout,*)' *** Warning ***'
146      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
147      WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
148      WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
149      abort_message='voir ci-dessus'
150      CALL abort_physic(modname,abort_message,1)
151    ENDIF
152
153    IF ( is_lic > is_sic ) THEN
154      WRITE(lunout,*)' *** Warning ***'
155      WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
156      WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
157      WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
158      abort_message='voir ci-dessus'
159      CALL abort_physic(modname,abort_message,1)
160    ENDIF
161
162!****************************************************************************************
163! Validation of ocean mode
164!
165!****************************************************************************************
166
167    IF (type_ocean /= 'slab  ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
168       WRITE(lunout,*)' *** Warning ***'
169       WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
170       abort_message='option pour l''ocean non valable'
171       CALL abort_physic(modname,abort_message,1)
172    ENDIF
173
174    iflag_pbl_surface_t2m_bug=0
175    CALL getin_p('iflag_pbl_surface_t2m_bug',iflag_pbl_surface_t2m_bug)
176    WRITE(lunout,*) 'iflag_pbl_surface_t2m_bug=',iflag_pbl_surface_t2m_bug
177!FC
178!    iflag_frein = 0
179!    CALL getin_p('iflag_frein',iflag_frein)
180!
181!jyg<
182!****************************************************************************************
183! Allocate variables for pbl splitting
184!
185!****************************************************************************************
186
187    CALL wx_pbl_init
188!>jyg
189
190  END SUBROUTINE pbl_surface_init
191
192#ifdef ISO
193  SUBROUTINE pbl_surface_init_iso(xtsnow_rst,Rland_ice_rst)
194
195! This routine should be called after the restart file has been read.
196! This routine initialize the restart variables and does some validation tests
197! for the index of the different surfaces and tests the choice of type of ocean.
198
199    USE indice_sol_mod
200    USE print_control_mod, ONLY: lunout
201#ifdef ISOVERIF
202    USE isotopes_mod, ONLY: iso_eau,ridicule
203    USE isotopes_verif_mod
204#endif
205    USE dimsoil_mod_h, ONLY: nsoilmx
206    IMPLICIT NONE
207 
208! Input variables
209!****************************************************************************************
210    REAL, DIMENSION(niso,klon, nbsrf), INTENT(IN)          :: xtsnow_rst
211    REAL, DIMENSION(niso,klon), INTENT(IN)          :: Rland_ice_rst
212 
213! Local variables
214!****************************************************************************************
215    INTEGER                       :: ierr
216    CHARACTER(len=80)             :: abort_message
217    CHARACTER(len = 20)           :: modname = 'pbl_surface_init'
218    integer i,ixt
219   
220!****************************************************************************************
221! Allocate and initialize module variables with fields read from restart file.
222!
223!****************************************************************************************   
224
225    ALLOCATE(xtsnow(niso,klon,nbsrf), stat=ierr)
226    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
227
228    ALLOCATE(Rland_ice(niso,klon), stat=ierr)
229    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
230
231    ALLOCATE(Roce(niso,klon), stat=ierr)
232    IF (ierr /= 0) CALL abort_physic('pbl_surface_init', 'pb in allocation',1)
233
234    xtsnow(:,:,:)  = xtsnow_rst(:,:,:)
235    Rland_ice(:,:) = Rland_ice_rst(:,:)
236    Roce(:,:)      = 0.0
237
238#ifdef ISOVERIF
239      IF (iso_eau >= 0) THEN
240         CALL iso_verif_egalite_vect2D( &
241     &           xtsnow,snow, &
242     &           'pbl_surface_mod 170',niso,klon,nbsrf)
243         DO i=1,klon 
244            IF (iso_eau >= 0) THEN 
245              CALL iso_verif_egalite(Rland_ice(iso_eau,i),1.0, &
246     &         'pbl_surf_mod 177')
247            ENDIF
248         ENDDO
249      ENDIF
250#endif
251
252  END SUBROUTINE pbl_surface_init_iso
253#endif
254
255
256!****************************************************************************************
257
258
259  SUBROUTINE pbl_surface( &
260       dtime,     date0,     itap,     jour,          &
261       debut,     lafin,                              &
262       rlon,      rlat,      rugoro,   rmu0,          &
263       lwdown_m,  cldt,                               &
264       rain_f,    snow_f,    bs_f, solsw_m,  solswfdiff_m, sollw_m,       &
265       gustiness,                                     &
266       t,         q,        qbs,  u,        v,        &
267!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
268!!       t_x,       q_x,       t_w,      q_w,           &
269       wake_dlt,             wake_dlq,                &
270       wake_cstar,           wake_s,                  &
271!!!
272       pplay,     paprs,     pctsrf,                  &
273       ts,SFRWL,   alb_dir, alb_dif,ustar, u10m, v10m,wstar, &
274       cdragh,    cdragm,   zu1,    zv1,              &
275!jyg<   (26/09/2019)
276       beta, &
277!>jyg
278       alb_dir_m,    alb_dif_m,  zxsens,   zxevap,  zxsnowerosion,      &
279       alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
280       zxtsol,    zxfluxlat, zt2m,     qsat2m, zn2mout,                 &
281       d_t,       d_q,    d_qbs,    d_u,      d_v, d_t_diss,            &
282!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
283       d_t_w,     d_q_w,                             &
284       d_t_x,     d_q_x,                             &
285!!       d_wake_dlt,d_wake_dlq,                         &
286       zxsens_x,  zxfluxlat_x,zxsens_w,zxfluxlat_w,  &
287!!!
288!!! nrlmd le 13/06/2011
289       delta_tsurf,wake_dens,cdragh_x,cdragh_w,      &
290       cdragm_x,cdragm_w,kh,kh_x,kh_w,               &
291!!!
292       zcoefh,    zcoefm,    slab_wfbils,            &
293       qsol,    zq2m,      s_pblh,   s_plcl,         &
294!!!
295!!! jyg le 08/02/2012
296       s_pblh_x, s_plcl_x,   s_pblh_w, s_plcl_w,     &
297!!!
298       s_capCL,   s_oliqCL,  s_cteiCL, s_pblT,       &
299       s_therm,   s_trmb1,   s_trmb2,  s_trmb3,      &
300       zustar,zu10m,  zv10m,    fder_print,          &
301       zxqsurf, delta_qsurf,                         &
302       rh2m,      zxfluxu,  zxfluxv,                 &
303       z0m, z0h,   agesno,  sollw,    solsw,         &
304       d_ts,      evap,    fluxlat,   t2m,           &
305       wfbils,    wfevap,                            &
306       flux_t,   flux_u, flux_v,                     &
307       dflux_t,   dflux_q,   zxsnow,                 &
308!jyg<
309!!       zxfluxt,   zxfluxq,   q2m,      flux_q, tke,   &
310       zxfluxt,   zxfluxq, zxfluxqbs,   q2m, flux_q, flux_qbs, tke_x, eps_x, &
311!>jyg
312!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
313!!        tke_x,     tke_w                              &
314       wake_dltke,                                   &
315        treedrg                                      &
316!FC
317!!!
318#ifdef ISO
319     &   ,xtrain_f, xtsnow_f,xt, &
320     &   wake_dlxt,zxxtevap,xtevap, &
321     &   d_xt,d_xt_w,d_xt_x, &
322     &   xtsol,dflux_xt,zxxtsnow,zxfluxxt,flux_xt, &
323     &   h1_diag,runoff_diag,xtrunoff_diag &
324#endif     
325     &   )
326!****************************************************************************************
327! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
328! Objet: interface de "couche limite" (diffusion verticale)
329!
330!AA REM:
331!AA-----
332!AA Tout ce qui a trait au traceurs est dans phytrac maintenant
333!AA pour l'instant le calcul de la couche limite pour les traceurs
334!AA se fait avec cltrac et ne tient pas compte de la differentiation
335!AA des sous-fraction de sol.
336!AA REM bis :
337!AA----------
338!AA Pour pouvoir extraire les coefficient d'echanges et le vent
339!AA dans la premiere couche, 3 champs supplementaires ont ete crees
340!AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
341!AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
342!AA si les informations des subsurfaces doivent etre prises en compte
343!AA il faudra sortir ces memes champs en leur ajoutant une dimension,
344!AA c'est a dire nbsrf (nbre de subsurface).
345!
346! Arguments:
347!
348! dtime----input-R- interval du temps (secondes)
349! itap-----input-I- numero du pas de temps
350! date0----input-R- jour initial
351! t--------input-R- temperature (K)
352! q--------input-R- vapeur d'eau (kg/kg)
353! u--------input-R- vitesse u
354! v--------input-R- vitesse v
355! wake_dlt-input-R- temperatre difference between (w) and (x) (K)
356! wake_dlq-input-R- humidity difference between (w) and (x) (kg/kg)
357!wake_cstar-input-R- wake gust front speed (m/s)
358! wake_s---input-R- wake fractionnal area
359! ts-------input-R- temperature du sol (en Kelvin)
360! paprs----input-R- pression a intercouche (Pa)
361! pplay----input-R- pression au milieu de couche (Pa)
362! rlat-----input-R- latitude en degree
363! z0m, z0h ----input-R- longeur de rugosite (en m)
364! Martin
365! cldt-----input-R- total cloud fraction
366! Martin
367!
368! d_t------output-R- le changement pour "t"
369! d_q------output-R- le changement pour "q"
370! d_u------output-R- le changement pour "u"
371! d_v------output-R- le changement pour "v"
372! d_ts-----output-R- le changement pour "ts"
373! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
374!                    (orientation positive vers le bas)
375! tke_x---input/output-R- tke in the (x) region (kg/m**2/s)
376! wake_dltke-input/output-R- tke difference between (w) and (x) (kg/m**2/s)
377! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
378! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
379! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
380! dflux_t--output-R- derive du flux sensible
381! dflux_q--output-R- derive du flux latent
382! zu1------output-R- le vent dans la premiere couche
383! zv1------output-R- le vent dans la premiere couche
384! trmb1----output-R- deep_cape
385! trmb2----output-R- inhibition
386! trmb3----output-R- Point Omega
387! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
388! plcl-----output-R- Niveau de condensation
389! pblh-----output-R- HCL
390! pblT-----output-R- T au nveau HCL
391! treedrg--output-R- tree drag (m)               
392!
393    USE carbon_cycle_mod,   ONLY : carbon_cycle_cpl, carbon_cycle_tr, level_coupling_esm
394    USE carbon_cycle_mod,   ONLY : co2_send, nbcf_out, fields_out, yfields_out, cfname_out
395    use hbtm_mod, only: hbtm
396    USE indice_sol_mod
397    USE time_phylmdz_mod,   ONLY : day_ini,annee_ref,itau_phy
398    USE mod_grid_phy_lmdz,  ONLY : nbp_lon, nbp_lat, grid1dto2d_glo
399    USE print_control_mod,  ONLY : prt_level,lunout
400#ifdef ISO
401  USE isotopes_mod, ONLY: Rdefault,iso_eau
402#ifdef ISOVERIF
403        USE isotopes_verif_mod
404#endif
405#ifdef ISOTRAC
406        USE isotrac_mod, only: index_iso
407#endif
408#endif
409USE dimpft_mod_h
410    USE flux_arp_mod_h
411    USE compbl_mod_h
412    USE yoethf_mod_h
413        USE clesphys_mod_h
414    USE ioipsl_getin_p_mod, ONLY : getin_p
415    use phys_state_var_mod, only: ds_ns, dt_ns, delta_sst, delta_sal, dter, &
416         dser, dt_ds, zsig, zmea
417    use phys_output_var_mod, only: tkt, tks, taur, sss
418    use lmdz_blowing_snow_ini, only : zeta_bs
419    use wxios_mod, ONLY: missing_val_xios => missing_val, using_xios
420    USE netcdf, only: missing_val_netcdf => nf90_fill_real
421    USE dimsoil_mod_h, ONLY: nsoilmx
422
423    USE yomcst_mod_h
424IMPLICIT NONE
425
426    INCLUDE "FCTTRE.h"
427!FC
428
429!****************************************************************************************
430    REAL,                         INTENT(IN)        :: dtime   ! time interval (s)
431    REAL,                         INTENT(IN)        :: date0   ! initial day
432    INTEGER,                      INTENT(IN)        :: itap    ! time step
433    INTEGER,                      INTENT(IN)        :: jour    ! current day of the year
434    LOGICAL,                      INTENT(IN)        :: debut   ! true if first run step
435    LOGICAL,                      INTENT(IN)        :: lafin   ! true if last run step
436    REAL, DIMENSION(klon),        INTENT(IN)        :: rlon    ! longitudes in degrees
437    REAL, DIMENSION(klon),        INTENT(IN)        :: rlat    ! latitudes in degrees
438    REAL, DIMENSION(klon),        INTENT(IN)        :: rugoro  ! rugosity length
439    REAL, DIMENSION(klon),        INTENT(IN)        :: rmu0    ! cosine of solar zenith angle
440    REAL, DIMENSION(klon),        INTENT(IN)        :: rain_f  ! rain fall
441    REAL, DIMENSION(klon),        INTENT(IN)        :: snow_f  ! snow fall
442    REAL, DIMENSION(klon),        INTENT(IN)        :: bs_f  ! blowing snow fall
443    REAL, DIMENSION(klon),        INTENT(IN)        :: solsw_m ! net shortwave radiation at mean surface
444    REAL, DIMENSION(klon),        INTENT(IN)        :: solswfdiff_m ! diffuse fraction fordownward shortwave radiation at mean surface
445    REAL, DIMENSION(klon),        INTENT(IN)        :: sollw_m ! net longwave radiation at mean surface
446    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t       ! temperature (K)
447    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q       ! water vapour (kg/kg)
448    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: qbs       ! blowing snow specific content (kg/kg)
449    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: u       ! u speed
450    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: v       ! v speed
451    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pplay   ! mid-layer pression (Pa)
452    REAL, DIMENSION(klon,klev+1), INTENT(IN)        :: paprs   ! pression between layers (Pa)
453    REAL, DIMENSION(klon, nbsrf), INTENT(IN)        :: pctsrf  ! sub-surface fraction
454! Martin
455    REAL, DIMENSION(klon),        INTENT(IN)        :: lwdown_m ! downward longwave radiation at mean s   
456    REAL, DIMENSION(klon),        INTENT(IN)        :: gustiness ! gustiness
457
458    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud
459
460#ifdef ISO
461    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: xt       ! water vapour (kg/kg)
462    REAL, DIMENSION(ntraciso,klon),        INTENT(IN)        :: xtrain_f  ! rain fall
463    REAL, DIMENSION(ntraciso,klon),        INTENT(IN)        :: xtsnow_f  ! snow fall
464#endif
465
466!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
467!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t_x       ! Temp\'erature hors poche froide
468!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: t_w       ! Temp\'erature dans la poches froide
469!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q_x       !
470!!    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: q_w       ! Pareil pour l'humidit\'e
471    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: wake_dlt  !temperature difference between (w) and (x) (K)
472    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: wake_dlq  !humidity difference between (w) and (x) (K)
473    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_s    ! Fraction de poches froides
474    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_cstar! Vitesse d'expansion des poches froides
475    REAL, DIMENSION(klon),        INTENT(IN)        :: wake_dens
476!!!
477#ifdef ISO
478    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)        :: wake_dlxt   
479#endif
480! Input/Output variables
481!****************************************************************************************
482!jyg<
483    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: beta    ! Aridity factor
484!>jyg
485    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ts      ! temperature at surface (K)
486    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: delta_tsurf !surface temperature difference between
487                                                                   !wake and off-wake regions
488!albedo SB >>>
489    REAL, DIMENSIOn(6),intent(in) :: SFRWL
490    REAL, DIMENSION(klon, nsw, nbsrf), INTENT(INOUT)     :: alb_dir,alb_dif
491!albedo SB <<<
492!jyg Pourquoi ustar et wstar sont-elles INOUT ?
493    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: ustar   ! u* (m/s)
494    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)   :: wstar   ! w* (m/s)
495    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: u10m    ! u speed at 10m
496    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: v10m    ! v speed at 10m
497!jyg<
498!!    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
499    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke_x
500!>jyg
501
502!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
503    REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: wake_dltke ! TKE_w - TKE_x
504!!!
505
506! Output variables
507!****************************************************************************************
508    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(OUT)   :: eps_x      ! TKE dissipation rate
509
510    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh     ! drag coefficient for T and Q
511    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm     ! drag coefficient for wind
512    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu1        ! u wind speed in first layer
513    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv1        ! v wind speed in first layer
514!albedo SB >>>
515    REAL, DIMENSION(klon, nsw),   INTENT(OUT)       :: alb_dir_m,alb_dif_m
516!albedo SB <<<
517    ! Martin
518    REAL, DIMENSION(klon),        INTENT(OUT)       :: alb3_lic
519    ! Martin
520    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens     ! sensible heat flux at surface with inversed sign
521                                                                  ! (=> positive sign upwards)
522    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxevap     ! water vapour flux at surface, positiv upwards
523    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsnowerosion     ! blowing snow flux at surface
524    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxtsol     ! temperature at surface, mean for each grid point
525!!! jyg le ???
526    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_t_w      !   !
527    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_q_w      !      !  Tendances dans les poches
528    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_t_x      !   !
529    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_q_x      !      !  Tendances hors des poches
530!!! jyg
531    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat  ! latent flux, mean for each grid point
532    REAL, DIMENSION(klon),        INTENT(OUT)       :: zt2m       ! temperature at 2m, mean for each grid point
533    INTEGER, DIMENSION(klon, 6),  INTENT(OUT)       :: zn2mout    ! number of times the 2m temperature is out of the [tsol,temp]
534    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsat2m
535    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t        ! change in temperature
536    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_t_diss       ! change in temperature
537    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_q        ! change in water vapour
538    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_u        ! change in u speed
539    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_v        ! change in v speed
540    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: d_qbs        ! change in blowing snow specific content
541
542
543    REAL, INTENT(OUT):: zcoefh(:, :, :) ! (klon, klev, nbsrf + 1)
544    ! coef for turbulent diffusion of T and Q, mean for each grid point
545
546    REAL, INTENT(OUT):: zcoefm(:, :, :) ! (klon, klev, nbsrf + 1)
547    ! coef for turbulent diffusion of U and V (?), mean for each grid point
548#ifdef ISO
549    REAL, DIMENSION(ntraciso,klon),        INTENT(OUT)       :: zxxtevap     ! water vapour flux at surface, positiv upwards
550    REAL, DIMENSION(ntraciso,klon, klev),  INTENT(OUT)       :: d_xt        ! change in water vapour
551    REAL, DIMENSION(klon),                 INTENT(OUT)       :: runoff_diag
552    REAL, DIMENSION(niso,klon),            INTENT(OUT)       :: xtrunoff_diag
553    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)       :: d_xt_w
554    REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)       :: d_xt_x
555#endif
556
557
558
559!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
560    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens_x   ! Flux sensible hors poche
561    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxsens_w   ! Flux sensible dans la poche
562    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat_x! Flux latent hors poche
563    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxfluxlat_w! Flux latent dans la poche
564!!    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_wake_dlt
565!!    REAL, DIMENSION(klon,klev),   INTENT(OUT)       :: d_wake_dlq
566
567! Output only for diagnostics
568    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh_x
569    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh_w
570    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm_x
571    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm_w
572    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh
573    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh_x
574    REAL, DIMENSION(klon),        INTENT(OUT)       :: kh_w
575!!!
576    REAL, DIMENSION(klon),        INTENT(OUT)       :: slab_wfbils! heat balance at surface only for slab at ocean points
577    REAL, DIMENSION(klon),        INTENT(OUT)       :: qsol     ! water height in the soil (mm)
578    REAL, DIMENSION(klon),        INTENT(OUT)       :: zq2m       ! water vapour at 2m, mean for each grid point
579    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh     ! height of the planetary boundary layer(HPBL)
580!!! jyg le 08/02/2012
581    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh_x   ! height of the PBL in the off-wake region
582    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblh_w   ! height of the PBL in the wake region
583!!!
584    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl     ! condensation level
585!!! jyg le 08/02/2012
586    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl_x   ! condensation level in the off-wake region
587    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_plcl_w   ! condensation level in the wake region
588!!!
589    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_capCL    ! CAPE of PBL
590    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_oliqCL   ! liquid water intergral of PBL
591    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_cteiCL   ! cloud top instab. crit. of PBL
592    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_pblT     ! temperature at PBLH
593    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_therm    ! thermal virtual temperature excess
594    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb1    ! deep cape, mean for each grid point
595    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb2    ! inhibition, mean for each grid point
596    REAL, DIMENSION(klon),        INTENT(OUT)       :: s_trmb3    ! point Omega, mean for each grid point
597    REAL, DIMENSION(klon),        INTENT(OUT)       :: zustar     ! u*
598    REAL, DIMENSION(klon),        INTENT(OUT)       :: zu10m      ! u speed at 10m, mean for each grid point
599    REAL, DIMENSION(klon),        INTENT(OUT)       :: zv10m      ! v speed at 10m, mean for each grid point
600    REAL, DIMENSION(klon),        INTENT(OUT)       :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
601    REAL, DIMENSION(klon),        INTENT(OUT)       :: zxqsurf    ! humidity at surface, mean for each grid point
602    REAL, DIMENSION(klon),        INTENT(OUT)       :: delta_qsurf! humidity difference at surface, mean for each grid point
603    REAL, DIMENSION(klon),        INTENT(OUT)       :: rh2m       ! relative humidity at 2m
604    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxu    ! u wind tension, mean for each grid point
605    REAL, DIMENSION(klon, klev),  INTENT(OUT)       :: zxfluxv    ! v wind tension, mean for each grid point
606    REAL, DIMENSION(klon, nbsrf+1), INTENT(INOUT)   :: z0m,z0h      ! rugosity length (m)
607    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: agesno   ! age of snow at surface
608    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: solsw      ! net shortwave radiation at surface
609    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: sollw      ! net longwave radiation at surface
610    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: d_ts       ! change in temperature at surface
611    REAL, DIMENSION(klon, nbsrf), INTENT(INOUT)     :: evap       ! evaporation at surface
612    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: fluxlat    ! latent flux
613    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: t2m        ! temperature at 2 meter height
614    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfbils     ! heat balance at surface
615    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)       :: wfevap     ! water balance (evap) at surface weighted by srf
616    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t     ! sensible heat flux (CpT) J/m**2/s (W/m**2)
617                                                                  ! positve orientation downwards
618    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u     ! u wind tension (kg m/s)/(m**2 s) or Pascal
619    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v     ! v wind tension (kg m/s)/(m**2 s) or Pascal
620!FC
621    REAL, DIMENSION(klon, klev, nbsrf), INTENT(INOUT) :: treedrg  ! tree drag (m)               
622#ifdef ISO       
623    REAL, DIMENSION(niso,klon),   INTENT(OUT)       :: xtsol      ! water height in the soil (mm)
624    REAL, DIMENSION(ntraciso,klon, nbsrf)           :: xtevap     ! evaporation at surface
625    REAL, DIMENSION(klon),        INTENT(OUT)       :: h1_diag    ! just diagnostic, not useful
626#endif
627
628
629! Output not needed
630    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_t    ! change of sensible heat flux
631    REAL, DIMENSION(klon),       INTENT(OUT)        :: dflux_q    ! change of water vapour flux
632    REAL, DIMENSION(klon),       INTENT(OUT)        :: zxsnow     ! snow at surface, mean for each grid point
633    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxt    ! sensible heat flux, mean for each grid point
634    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxq    ! water vapour flux, mean for each grid point
635    REAL, DIMENSION(klon, klev), INTENT(OUT)        :: zxfluxqbs    ! blowing snow flux, mean for each grid point
636    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m        ! water vapour at 2 meter height
637    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q     ! water vapour flux(latent flux) (kg/m**2/s)
638    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_qbs   ! blowind snow vertical flux (kg/m**2
639
640#ifdef ISO   
641    REAL, DIMENSION(ntraciso,klon),              INTENT(OUT) :: dflux_xt    ! change of water vapour flux
642    REAL, DIMENSION(niso,klon),                  INTENT(OUT) :: zxxtsnow    ! snow at surface, mean for each grid point
643    REAL, DIMENSION(ntraciso,klon, klev),        INTENT(OUT) :: zxfluxxt    ! water vapour flux, mean for each grid point
644    REAL, DIMENSION(ntraciso,klon, klev, nbsrf), INTENT(OUT) :: flux_xt     ! water vapour flux(latent flux) (kg/m**2/s) 
645#endif
646
647! Martin
648! inlandsis
649    REAL, DIMENSION(klon),       INTENT(OUT)        :: qsnow      ! snow water content
650    REAL, DIMENSION(klon),       INTENT(OUT)        :: snowhgt    ! snow height
651    REAL, DIMENSION(klon),       INTENT(OUT)        :: to_ice     ! snow passed to ice
652    REAL, DIMENSION(klon),       INTENT(OUT)        :: sissnow    ! snow in snow model
653    REAL, DIMENSION(klon),       INTENT(OUT)        :: runoff     ! runoff on land ice
654! Martin
655
656! Local variables with attribute SAVE
657!****************************************************************************************
658    INTEGER, SAVE                            :: nhoridbg, nidbg   ! variables for IOIPSL
659!$OMP THREADPRIVATE(nhoridbg, nidbg)
660    LOGICAL, SAVE                            :: debugindex=.FALSE.
661!$OMP THREADPRIVATE(debugindex)
662    LOGICAL, SAVE                            :: first_call=.TRUE.
663!$OMP THREADPRIVATE(first_call)
664    CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
665!$OMP THREADPRIVATE(cl_surf)
666    REAL, SAVE                               :: beta_land         ! beta for wx_dts
667!$OMP THREADPRIVATE(beta_land)
668
669! Other local variables
670!****************************************************************************************
671! >> PC
672    INTEGER                            :: ierr
673    INTEGER                            :: n
674! << PC
675    INTEGER                            :: iflag_split, iflag_split_ref
676    INTEGER                            :: i, k, nsrf
677    INTEGER                            :: knon, j
678    INTEGER                            :: idayref
679    INTEGER , DIMENSION(klon)          :: ni
680    REAL                               :: yt1_new
681    REAL                               :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
682    REAL                               :: amn, amx
683    REAL                               :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
684    REAL, DIMENSION(klon)              :: r_co2_ppm     ! taux CO2 atmosphere
685    REAL, DIMENSION(klon)              :: yts, yz0m, yz0h, ypct
686    REAL, DIMENSION(klon)              :: yz0h_old
687!albedo SB >>>
688    REAL, DIMENSION(klon)              :: yalb,yalb_vis
689!albedo SB <<<
690    REAL, DIMENSION(klon)              :: yt1, yq1, yu1, yv1, yqbs1
691    REAL, DIMENSION(klon)              :: yqa
692    REAL, DIMENSION(klon)              :: ysnow, yqsurf, yagesno, yqsol
693    REAL, DIMENSION(klon)              :: yrain_f, ysnow_f, ybs_f
694#ifdef ISO
695    REAL, DIMENSION(ntraciso,klon)     :: yxt1
696    REAL, DIMENSION(niso,klon)         :: yxtsnow, yxtsol   
697    REAL, DIMENSION(ntraciso,klon)     :: yxtrain_f, yxtsnow_f
698    REAL, DIMENSION(klon)              :: yrunoff_diag
699    REAL, DIMENSION(niso,klon)         :: yxtrunoff_diag
700    REAL, DIMENSION(niso,klon)         :: yRland_ice   
701#endif
702    REAL, DIMENSION(klon)              :: ysolsw, ysollw
703    REAL, DIMENSION(klon)              :: yfder
704    REAL, DIMENSION(klon)              :: yrugoro
705    REAL, DIMENSION(klon)              :: yfluxlat
706    REAL, DIMENSION(klon)              :: yfluxbs
707    REAL, DIMENSION(klon)              :: y_d_ts
708    REAL, DIMENSION(klon)              :: y_flux_t1, y_flux_q1
709    REAL, DIMENSION(klon)              :: y_dflux_t, y_dflux_q
710#ifdef ISO
711    REAL, DIMENSION(ntraciso,klon)     ::  y_flux_xt1
712    REAL, DIMENSION(ntraciso,klon)     ::  y_dflux_xt
713#endif
714    REAL, DIMENSION(klon)              :: y_flux_u1, y_flux_v1
715    REAL, DIMENSION(klon)              :: y_flux_bs, y_flux0
716    REAL, DIMENSION(klon)              :: yt2m, yq2m, yu10m
717    INTEGER, DIMENSION(klon, nbsrf, 6) :: yn2mout, yn2mout_x, yn2mout_w
718    INTEGER, DIMENSION(klon, nbsrf, 6) :: n2mout, n2mout_x, n2mout_w
719    REAL, DIMENSION(klon)              :: yustar
720    REAL, DIMENSION(klon)              :: ywstar
721    REAL, DIMENSION(klon)              :: ywindsp
722    REAL, DIMENSION(klon)              :: yt10m, yq10m
723    REAL, DIMENSION(klon)              :: ypblh
724    REAL, DIMENSION(klon)              :: ylcl
725    REAL, DIMENSION(klon)              :: ycapCL
726    REAL, DIMENSION(klon)              :: yoliqCL
727    REAL, DIMENSION(klon)              :: ycteiCL
728    REAL, DIMENSION(klon)              :: ypblT
729    REAL, DIMENSION(klon)              :: ytherm
730    REAL, DIMENSION(klon)              :: ytrmb1
731    REAL, DIMENSION(klon)              :: ytrmb2
732    REAL, DIMENSION(klon)              :: ytrmb3
733    REAL, DIMENSION(klon)              :: uzon, vmer
734    REAL, DIMENSION(klon)              :: tair1, qair1, tairsol
735    REAL, DIMENSION(klon)              :: psfce, patm
736    REAL, DIMENSION(klon)              :: qairsol, zgeo1, speed, zri1, pref !speed, zri1, pref, added by Fuxing WANG, 04/03/2015
737    REAL, DIMENSION(klon)              :: yz0h_oupas
738    REAL, DIMENSION(klon)              :: yfluxsens
739    REAL, DIMENSION(klon)              :: AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0
740    REAL, DIMENSION(klon)              :: AcoefH, AcoefQ, BcoefH, BcoefQ
741#ifdef ISO
742    REAL, DIMENSION(ntraciso,klon)     :: AcoefXT, BcoefXT
743#endif
744    REAL, DIMENSION(klon)              :: AcoefU, AcoefV, BcoefU, BcoefV
745    REAL, DIMENSION(klon)              :: AcoefQBS, BcoefQBS
746    REAL, DIMENSION(klon)              :: ypsref
747    REAL, DIMENSION(klon)              :: yevap, yevap_pot, ytsurf_new, yalb3_new
748!albedo SB >>>
749    REAL, DIMENSION(klon,nsw)          :: yalb_dir_new, yalb_dif_new
750!albedo SB <<<
751    REAL, DIMENSION(klon)              :: ztsol
752    REAL, DIMENSION(klon)              :: meansqT ! mean square deviation of subsurface temperatures
753    REAL, DIMENSION(klon)              :: alb_m  ! mean albedo for whole SW interval
754    REAL, DIMENSION(klon,klev)         :: y_d_t, y_d_q, y_d_t_diss, y_d_qbs
755    REAL, DIMENSION(klon,klev)         :: y_d_u, y_d_v
756    REAL, DIMENSION(klon,klev)         :: y_flux_t, y_flux_q, y_flux_qbs
757    REAL, DIMENSION(klon,klev)         :: y_flux_u, y_flux_v
758    REAL, DIMENSION(klon,klev)         :: ycoefh,ycoefm,ycoefq,ycoefqbs
759    REAL, DIMENSION(klon)              :: ycdragh, ycdragq, ycdragm
760    REAL, DIMENSION(klon,klev)         :: yu, yv
761    REAL, DIMENSION(klon,klev)         :: yt, yq, yqbs
762#ifdef ISO
763    REAL, DIMENSION(ntraciso,klon)      :: yxtevap
764    REAL, DIMENSION(ntraciso,klon,klev) :: y_d_xt
765    REAL, DIMENSION(ntraciso,klon,klev) :: y_flux_xt
766    REAL, DIMENSION(ntraciso,klon,klev) :: yxt   
767#endif
768    REAL, DIMENSION(klon,klev)         :: ypplay, ydelp
769    REAL, DIMENSION(klon,klev)         :: delp
770    REAL, DIMENSION(klon,klev+1)       :: ypaprs
771    REAL, DIMENSION(klon,klev+1)       :: ytke, yeps
772    REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
773!FC
774    REAL, DIMENSION(klon,nvm_lmdz)          :: yveget
775    REAL, DIMENSION(klon,nvm_lmdz)          :: ylai
776    REAL, DIMENSION(klon,nvm_lmdz)          :: yheight
777    REAL, DIMENSION(klon,klev)              :: y_d_u_frein
778    REAL, DIMENSION(klon,klev)              :: y_d_v_frein
779    REAL, DIMENSION(klon,klev)              :: y_treedrg
780!FC
781
782    CHARACTER(len=80)                  :: abort_message
783    CHARACTER(len=20)                  :: modname = 'pbl_surface'
784    LOGICAL, PARAMETER                 :: zxli=.FALSE. ! utiliser un jeu de fonctions simples
785    LOGICAL, PARAMETER                 :: check=.FALSE.
786
787!!! nrlmd le 02/05/2011
788!!! jyg le 07/02/2012
789    REAL, DIMENSION(klon)              :: ywake_s, ywake_cstar, ywake_dens
790!!!
791    REAL, DIMENSION(klon,klev+1)       :: ytke_x, ytke_w, yeps_x, yeps_w
792    REAL, DIMENSION(klon,klev+1)       :: ywake_dltke
793    REAL, DIMENSION(klon,klev)         :: yu_x, yv_x, yu_w, yv_w
794    REAL, DIMENSION(klon,klev)         :: yt_x, yq_x, yt_w, yq_w
795    REAL, DIMENSION(klon,klev)         :: ycoefh_x, ycoefm_x, ycoefh_w, ycoefm_w
796    REAL, DIMENSION(klon,klev)         :: ycoefq_x, ycoefq_w
797    REAL, DIMENSION(klon)              :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
798    REAL, DIMENSION(klon)              :: ycdragm_x, ycdragm_w
799    REAL, DIMENSION(klon)              :: AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x
800    REAL, DIMENSION(klon)              :: AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w
801    REAL, DIMENSION(klon)              :: AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x
802    REAL, DIMENSION(klon)              :: AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w
803    REAL, DIMENSION(klon)              :: y_flux_t1_x, y_flux_q1_x, y_flux_t1_w, y_flux_q1_w
804    REAL, DIMENSION(klon)              :: y_flux_u1_x, y_flux_v1_x, y_flux_u1_w, y_flux_v1_w
805    REAL, DIMENSION(klon,klev)         :: y_flux_t_x, y_flux_q_x, y_flux_t_w, y_flux_q_w
806    REAL, DIMENSION(klon,klev)         :: y_flux_u_x, y_flux_v_x, y_flux_u_w, y_flux_v_w
807    REAL, DIMENSION(klon)              :: yfluxlat_x, yfluxlat_w
808    REAL, DIMENSION(klon,klev)         :: y_d_t_x, y_d_q_x, y_d_t_w, y_d_q_w
809    REAL, DIMENSION(klon,klev)         :: y_d_t_diss_x, y_d_t_diss_w
810    REAL, DIMENSION(klon,klev)         :: d_t_diss_x, d_t_diss_w
811    REAL, DIMENSION(klon,klev)         :: y_d_u_x, y_d_v_x, y_d_u_w, y_d_v_w
812    REAL, DIMENSION(klon, klev, nbsrf) :: flux_t_x, flux_q_x, flux_t_w, flux_q_w
813    REAL, DIMENSION(klon, klev, nbsrf) :: flux_u_x, flux_v_x, flux_u_w, flux_v_w
814    REAL, DIMENSION(klon, nbsrf)       :: fluxlat_x, fluxlat_w
815    REAL, DIMENSION(klon, klev)        :: zxfluxt_x, zxfluxq_x, zxfluxt_w, zxfluxq_w
816    REAL, DIMENSION(klon, klev)        :: zxfluxu_x, zxfluxv_x, zxfluxu_w, zxfluxv_w
817    REAL                               :: zx_qs_surf, zcor_surf, zdelta_surf
818!jyg<
819    REAL, DIMENSION(klon)              :: ybeta
820    REAL, DIMENSION(klon)              :: ybeta_prev
821!>jyg
822    REAL, DIMENSION(klon, klev)        :: d_u_x
823    REAL, DIMENSION(klon, klev)        :: d_u_w
824    REAL, DIMENSION(klon, klev)        :: d_v_x
825    REAL, DIMENSION(klon, klev)        :: d_v_w
826
827    REAL, DIMENSION(klon,klev)         :: CcoefH, CcoefQ, DcoefH, DcoefQ
828    REAL, DIMENSION(klon,klev)         :: CcoefU, CcoefV, DcoefU, DcoefV
829    REAL, DIMENSION(klon,klev)         :: CcoefQBS, DcoefQBS
830    REAL, DIMENSION(klon,klev)         :: CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x
831    REAL, DIMENSION(klon,klev)         :: CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w
832    REAL, DIMENSION(klon,klev)         :: CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x
833    REAL, DIMENSION(klon,klev)         :: CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w
834    REAL, DIMENSION(klon,klev)         :: Kcoef_hq, Kcoef_m, gama_h, gama_q
835    REAL, DIMENSION(klon,klev)         :: gama_qbs, Kcoef_qbs
836    REAL, DIMENSION(klon,klev)         :: Kcoef_hq_x, Kcoef_m_x, gama_h_x, gama_q_x
837    REAL, DIMENSION(klon,klev)         :: Kcoef_hq_w, Kcoef_m_w, gama_h_w, gama_q_w
838    REAL, DIMENSION(klon)              :: alf_1, alf_2, alf_1_x, alf_2_x, alf_1_w, alf_2_w
839#ifdef ISO
840    REAL, DIMENSION(ntraciso,klon,klev)         :: yxt_x, yxt_w
841    REAL, DIMENSION(ntraciso,klon)              :: y_flux_xt1_x , y_flux_xt1_w   
842    REAL, DIMENSION(ntraciso,klon,klev)         :: y_flux_xt_x,y_d_xt_x,zxfluxxt_x
843    REAL, DIMENSION(ntraciso,klon,klev)         :: y_flux_xt_w,y_d_xt_w,zxfluxxt_w
844    REAL, DIMENSION(ntraciso,klon,klev,nbsrf)   :: flux_xt_x, flux_xt_w
845    REAL, DIMENSION(ntraciso,klon)              :: AcoefXT_x, BcoefXT_x
846    REAL, DIMENSION(ntraciso,klon)              :: AcoefXT_w, BcoefXT_w
847    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT, DcoefXT
848    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT_x, DcoefXT_x
849    REAL, DIMENSION(ntraciso,klon,klev)         :: CcoefXT_w, DcoefXT_w
850    REAL, DIMENSION(ntraciso,klon,klev)         :: gama_xt,gama_xt_x,gama_xt_w
851#endif
852!!!
853!!!jyg le 08/02/2012
854    REAL, DIMENSION(klon, nbsrf)       :: windsp
855!
856    REAL, DIMENSION(klon, nbsrf)       :: t2m_x
857    REAL, DIMENSION(klon, nbsrf)       :: q2m_x
858    REAL, DIMENSION(klon)              :: rh2m_x
859    REAL, DIMENSION(klon)              :: qsat2m_x
860    REAL, DIMENSION(klon, nbsrf)       :: u10m_x
861    REAL, DIMENSION(klon, nbsrf)       :: v10m_x
862    REAL, DIMENSION(klon, nbsrf)       :: ustar_x
863    REAL, DIMENSION(klon, nbsrf)       :: wstar_x
864!             
865    REAL, DIMENSION(klon, nbsrf)       :: pblh_x
866    REAL, DIMENSION(klon, nbsrf)       :: plcl_x
867    REAL, DIMENSION(klon, nbsrf)       :: capCL_x
868    REAL, DIMENSION(klon, nbsrf)       :: oliqCL_x
869    REAL, DIMENSION(klon, nbsrf)       :: cteiCL_x
870    REAL, DIMENSION(klon, nbsrf)       :: pblt_x
871    REAL, DIMENSION(klon, nbsrf)       :: therm_x
872    REAL, DIMENSION(klon, nbsrf)       :: trmb1_x
873    REAL, DIMENSION(klon, nbsrf)       :: trmb2_x
874    REAL, DIMENSION(klon, nbsrf)       :: trmb3_x
875!
876    REAL, DIMENSION(klon, nbsrf)       :: t2m_w
877    REAL, DIMENSION(klon, nbsrf)       :: q2m_w
878    REAL, DIMENSION(klon)              :: rh2m_w
879    REAL, DIMENSION(klon)              :: qsat2m_w
880    REAL, DIMENSION(klon, nbsrf)       :: u10m_w
881    REAL, DIMENSION(klon, nbsrf)       :: v10m_w
882    REAL, DIMENSION(klon, nbsrf)       :: ustar_w
883    REAL, DIMENSION(klon, nbsrf)       :: wstar_w
884!                           
885    REAL, DIMENSION(klon, nbsrf)       :: pblh_w
886    REAL, DIMENSION(klon, nbsrf)       :: plcl_w
887    REAL, DIMENSION(klon, nbsrf)       :: capCL_w
888    REAL, DIMENSION(klon, nbsrf)       :: oliqCL_w
889    REAL, DIMENSION(klon, nbsrf)       :: cteiCL_w
890    REAL, DIMENSION(klon, nbsrf)       :: pblt_w
891    REAL, DIMENSION(klon, nbsrf)       :: therm_w
892    REAL, DIMENSION(klon, nbsrf)       :: trmb1_w
893    REAL, DIMENSION(klon, nbsrf)       :: trmb2_w
894    REAL, DIMENSION(klon, nbsrf)       :: trmb3_w
895!
896    REAL, DIMENSION(klon)       :: yt2m_x
897    REAL, DIMENSION(klon)       :: yq2m_x
898    REAL, DIMENSION(klon)       :: yt10m_x
899    REAL, DIMENSION(klon)       :: yq10m_x
900    REAL, DIMENSION(klon)       :: yu10m_x
901    REAL, DIMENSION(klon)       :: yv10m_x
902    REAL, DIMENSION(klon)       :: yustar_x
903    REAL, DIMENSION(klon)       :: ywstar_x
904!             
905    REAL, DIMENSION(klon)       :: ypblh_x
906    REAL, DIMENSION(klon)       :: ylcl_x
907    REAL, DIMENSION(klon)       :: ycapCL_x
908    REAL, DIMENSION(klon)       :: yoliqCL_x
909    REAL, DIMENSION(klon)       :: ycteiCL_x
910    REAL, DIMENSION(klon)       :: ypblt_x
911    REAL, DIMENSION(klon)       :: ytherm_x
912    REAL, DIMENSION(klon)       :: ytrmb1_x
913    REAL, DIMENSION(klon)       :: ytrmb2_x
914    REAL, DIMENSION(klon)       :: ytrmb3_x
915!
916    REAL, DIMENSION(klon)       :: yt2m_w
917    REAL, DIMENSION(klon)       :: yq2m_w
918    REAL, DIMENSION(klon)       :: yt10m_w
919    REAL, DIMENSION(klon)       :: yq10m_w
920    REAL, DIMENSION(klon)       :: yu10m_w
921    REAL, DIMENSION(klon)       :: yv10m_w
922    REAL, DIMENSION(klon)       :: yustar_w
923    REAL, DIMENSION(klon)       :: ywstar_w
924!                       
925    REAL, DIMENSION(klon)       :: ypblh_w
926    REAL, DIMENSION(klon)       :: ylcl_w
927    REAL, DIMENSION(klon)       :: ycapCL_w
928    REAL, DIMENSION(klon)       :: yoliqCL_w
929    REAL, DIMENSION(klon)       :: ycteiCL_w
930    REAL, DIMENSION(klon)       :: ypblt_w
931    REAL, DIMENSION(klon)       :: ytherm_w
932    REAL, DIMENSION(klon)       :: ytrmb1_w
933    REAL, DIMENSION(klon)       :: ytrmb2_w
934    REAL, DIMENSION(klon)       :: ytrmb3_w
935!
936    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
937    REAL, DIMENSION(klon)       :: zgeo1_x, tair1_x, qair1_x, tairsol_x
938!
939    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
940    REAL, DIMENSION(klon)       :: zgeo1_w, tair1_w, qair1_w, tairsol_w
941    REAL, DIMENSION(klon)       :: yus0, yvs0
942
943!!! jyg le 25/03/2013
944!!    Variables intermediaires pour le raccord des deux colonnes \`a la surface
945!jyg<
946!!    REAL   ::   dd_Ch
947!!    REAL   ::   dd_Cm
948!!    REAL   ::   dd_Kh
949!!    REAL   ::   dd_Km
950!!    REAL   ::   dd_u
951!!    REAL   ::   dd_v
952!!    REAL   ::   dd_t
953!!    REAL   ::   dd_q
954!!    REAL   ::   dd_AH
955!!    REAL   ::   dd_AQ
956!!    REAL   ::   dd_AU
957!!    REAL   ::   dd_AV
958!!    REAL   ::   dd_BH
959!!    REAL   ::   dd_BQ
960!!    REAL   ::   dd_BU
961!!    REAL   ::   dd_BV
962!!
963!!    REAL   ::   dd_KHp
964!!    REAL   ::   dd_KQp
965!!    REAL   ::   dd_KUp
966!!    REAL   ::   dd_KVp
967!>jyg
968
969!!!
970!!! nrlmd le 13/06/2011
971    REAL, DIMENSION(klon)              :: y_delta_flux_t1, y_delta_flux_q1, y_delta_flux_u1, y_delta_flux_v1
972    REAL, DIMENSION(klon)              :: y_delta_tsurf, y_delta_tsurf_new
973    REAL, DIMENSION(klon)              :: delta_coef, tau_eq
974    REAL, DIMENSION(klon)              :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
975    REAL, DIMENSION(klon)              :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
976    REAL, DIMENSION(klon)              :: y_delta_qsurf
977    REAL, DIMENSION(klon)              :: y_delta_qsats
978    REAL, DIMENSION(klon)              :: yg_T, yg_Q
979    REAL, DIMENSION(klon)              :: yGamma_dTs_phiT, yGamma_dQs_phiQ
980    REAL, DIMENSION(klon)              :: ydTs_ins, ydqs_ins
981!
982    REAL, PARAMETER                    :: facteur=2./sqrt(3.14)
983    REAL, PARAMETER                    :: inertia=2000.
984    REAL, DIMENSION(klon)              :: ydtsurf_th
985    REAL                               :: zdelta_surf_x,zdelta_surf_w,zx_qs_surf_x,zx_qs_surf_w
986    REAL                               :: zcor_surf_x,zcor_surf_w
987    REAL                               :: mod_wind_x, mod_wind_w
988    REAL                               :: rho1
989    REAL, DIMENSION(klon)              :: Kech_h           ! Coefficient d'echange pour l'energie
990    REAL, DIMENSION(klon)              :: Kech_h_x, Kech_h_w
991    REAL, DIMENSION(klon)              :: Kech_m
992    REAL, DIMENSION(klon)              :: Kech_m_x, Kech_m_w
993    REAL, DIMENSION(klon)              :: yts_x, yts_w
994    REAL, DIMENSION(klon)              :: yqsatsrf0_x, yqsatsrf0_w
995    REAL, DIMENSION(klon)              :: yqsurf_x, yqsurf_w
996!jyg<
997!!    REAL, DIMENSION(klon)              :: Kech_Hp, Kech_H_xp, Kech_H_wp
998!!    REAL, DIMENSION(klon)              :: Kech_Qp, Kech_Q_xp, Kech_Q_wp
999!!    REAL, DIMENSION(klon)              :: Kech_Up, Kech_U_xp, Kech_U_wp
1000!!    REAL, DIMENSION(klon)              :: Kech_Vp, Kech_V_xp, Kech_V_wp
1001!>jyg
1002
1003    REAL                               :: fact_cdrag
1004    REAL                               :: z1lay
1005
1006    REAL                               :: vent
1007!
1008! For debugging with IOIPSL
1009    INTEGER, DIMENSION(nbp_lon*nbp_lat)    :: ndexbg
1010    REAL                               :: zjulian
1011    REAL, DIMENSION(klon)              :: tabindx
1012    REAL, DIMENSION(nbp_lon,nbp_lat)         :: zx_lon, zx_lat
1013    REAL, DIMENSION(nbp_lon,nbp_lat)         :: debugtab
1014
1015
1016    REAL, DIMENSION(klon,nbsrf)        :: pblh         ! height of the planetary boundary layer
1017    REAL, DIMENSION(klon,nbsrf)        :: plcl         ! condensation level
1018    REAL, DIMENSION(klon,nbsrf)        :: capCL
1019    REAL, DIMENSION(klon,nbsrf)        :: oliqCL
1020    REAL, DIMENSION(klon,nbsrf)        :: cteiCL
1021    REAL, DIMENSION(klon,nbsrf)        :: pblT
1022    REAL, DIMENSION(klon,nbsrf)        :: therm
1023    REAL, DIMENSION(klon,nbsrf)        :: trmb1        ! deep cape
1024    REAL, DIMENSION(klon,nbsrf)        :: trmb2        ! inhibition
1025    REAL, DIMENSION(klon,nbsrf)        :: trmb3        ! point Omega
1026    REAL, DIMENSION(klon,nbsrf)        :: zx_rh2m, zx_qsat2m
1027    REAL, DIMENSION(klon,nbsrf)        :: zx_t1
1028    REAL, DIMENSION(klon, nbsrf)       :: alb          ! mean albedo for whole SW interval
1029    REAL, DIMENSION(klon,nbsrf)        :: snowerosion   
1030    REAL, DIMENSION(klon)              :: ylwdown      ! jg : temporary (ysollwdown)
1031    REAL, DIMENSION(klon)              :: ygustiness      ! jg : temporary (ysollwdown)
1032
1033    REAL                               :: zx_qs1, zcor1, zdelta1
1034
1035    ! Martin
1036    REAL, DIMENSION(klon, nbsrf)       :: sollwd ! net longwave radiation at surface
1037    REAL, DIMENSION(klon)              :: ytoice
1038    REAL, DIMENSION(klon)              :: ysnowhgt, yqsnow, ysissnow, yrunoff
1039    REAL, DIMENSION(klon)              :: yzmea
1040    REAL, DIMENSION(klon)              :: yzsig
1041    REAL, DIMENSION(klon)              :: ycldt
1042    REAL, DIMENSION(klon)              :: yrmu0
1043    ! Martin
1044    REAL, DIMENSION(klon)              :: yri0
1045
1046    REAL, DIMENSION(klon):: ydelta_sst, ydelta_sal, yds_ns, ydt_ns, ydter, &
1047         ydser, ydt_ds, ytkt, ytks, ytaur, ysss
1048    ! compression of delta_sst, delta_sal, ds_ns, dt_ns, dter, dser,
1049    ! dt_ds, tkt, tks, taur, sss on ocean points
1050    REAL :: missing_val
1051#ifdef ISO
1052    REAL, DIMENSION(klon)       :: h1
1053    INTEGER                     :: ixt
1054!#ifdef ISOVERIF
1055!    integer iso_verif_positif_nostop
1056!#endif   
1057#endif
1058
1059!****************************************************************************************
1060! End of declarations
1061!****************************************************************************************
1062      IF (using_xios) THEN
1063        missing_val=missing_val_xios
1064      ELSE
1065        missing_val=missing_val_netcdf
1066      ENDIF
1067
1068      IF (prt_level >=10) print *,' -> pbl_surface, itap ',itap
1069!
1070!!jyg      iflag_split = mod(iflag_pbl_split,2)
1071!!jyg      iflag_split = mod(iflag_pbl_split,10)
1072!
1073! Flags controlling the splitting of the turbulent boundary layer:
1074!   iflag_split_ref = 0  ==> no splitting
1075!                   = 1  ==> splitting without coupling with surface temperature
1076!                   = 2  ==> splitting with coupling with surface temperature over land
1077!                   = 3  ==> splitting over ocean; no splitting over land
1078!   iflag_split: actual flag controlling the splitting.
1079!   iflag_split = iflag_split_ref outside the sub-surface loop
1080!               = iflag_split_ref if iflag_split_ref = 0, 1, or 2
1081!               = 0 over land  if iflga_split_ref = 3
1082!               = 1 over ocean if iflga_split_ref = 3
1083
1084      iflag_split_ref = mod(iflag_pbl_split,10)
1085      iflag_split = iflag_split_ref
1086
1087#ifdef ISO     
1088#ifdef ISOVERIF
1089      DO i=1,klon
1090        DO ixt=1,niso
1091          CALL iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 608')
1092        ENDDO
1093      ENDDO
1094#endif
1095#ifdef ISOVERIF
1096      DO i=1,klon 
1097        IF (iso_eau >= 0) THEN 
1098          CALL iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
1099     &         'pbl_surf_mod 585',errmax,errmaxrel)
1100          CALL iso_verif_egalite_choix(xtsnow_f(iso_eau,i),snow_f(i), &
1101     &         'pbl_surf_mod 594',errmax,errmaxrel)
1102          IF (iso_verif_egalite_choix_nostop(xtsol(iso_eau,i),qsol(i), &
1103     &         'pbl_surf_mod 596',errmax,errmaxrel) == 1) THEN
1104                WRITE(*,*) 'i=',i
1105                STOP
1106          ENDIF
1107          DO nsrf=1,nbsrf
1108            CALL iso_verif_egalite_choix(xtsnow(iso_eau,i,nsrf),snow(i,nsrf), &
1109     &         'pbl_surf_mod 598',errmax,errmaxrel)
1110          ENDDO
1111        ENDIF !IF (iso_eau >= 0) THEN   
1112      ENDDO !DO i=1,knon 
1113      DO k=1,klev
1114        DO i=1,klon 
1115          IF (iso_eau >= 0) THEN 
1116            CALL iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &
1117     &           'pbl_surf_mod 595',errmax,errmaxrel)
1118          ENDIF !IF (iso_eau >= 0) THEN 
1119        ENDDO !DO i=1,knon 
1120      ENDDO !DO k=1,klev
1121#endif
1122#endif
1123
1124
1125!****************************************************************************************
1126! 1) Initialisation and validation tests
1127!    Only done first time entering this subroutine
1128!
1129!****************************************************************************************
1130
1131    IF (first_call) THEN
1132
1133       iflag_new_t2mq2m=1
1134       CALL getin_p('iflag_new_t2mq2m',iflag_new_t2mq2m)
1135       WRITE(lunout,*) 'pbl_iflag_new_t2mq2m=',iflag_new_t2mq2m
1136
1137       ok_bug_zg_wk_pbl=.TRUE.
1138       CALL getin_p('ok_bug_zg_wk_pbl',ok_bug_zg_wk_pbl)
1139       WRITE(lunout,*) 'ok_bug_zg_wk_pbl=',ok_bug_zg_wk_pbl
1140
1141       print*,'PBL SURFACE AVEC GUSTINESS'
1142       first_call=.FALSE.
1143     
1144       ! Initialize ok_flux_surf (for 1D model)
1145       IF (klon_glo>1) ok_flux_surf=.FALSE.
1146       IF (klon_glo>1) ok_forc_tsurf=.FALSE.
1147
1148       ! intialize beta_land
1149       beta_land = 0.5
1150       call getin_p('beta_land', beta_land)
1151       
1152       ! Initilize debug IO
1153       IF (debugindex .AND. mpi_size==1) THEN
1154          ! initialize IOIPSL output
1155          idayref = day_ini
1156          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
1157          CALL grid1dTo2d_glo(rlon,zx_lon)
1158          DO i = 1, nbp_lon
1159             zx_lon(i,1) = rlon(i+1)
1160             zx_lon(i,nbp_lat) = rlon(i+1)
1161          ENDDO
1162          CALL grid1dTo2d_glo(rlat,zx_lat)
1163          CALL histbeg("sous_index",nbp_lon,zx_lon(:,1),nbp_lat,zx_lat(1,:), &
1164               1,nbp_lon,1,nbp_lat, &
1165               itau_phy,zjulian,dtime,nhoridbg,nidbg)
1166          ! no vertical axis
1167          cl_surf(1)='ter'
1168          cl_surf(2)='lic'
1169          cl_surf(3)='oce'
1170          cl_surf(4)='sic'
1171          DO nsrf=1,nbsrf
1172             CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",nbp_lon, &
1173                  nbp_lat,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1174          ENDDO
1175
1176          CALL histend(nidbg)
1177          CALL histsync(nidbg)
1178
1179       ENDIF
1180       
1181    ENDIF
1182         
1183!****************************************************************************************
1184! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
1185! instead of ORCHIDEE)
1186    IF (qsol0>=0.) THEN
1187      PRINT*,'WARNING : On impose qsol=',qsol0
1188      qsol(:)=qsol0
1189#ifdef ISO
1190      DO ixt=1,niso
1191        xtsol(ixt,:)=qsol0*Rdefault(ixt)
1192      ENDDO
1193#ifdef ISOTRAC     
1194      DO ixt=1+niso,ntraciso
1195        xtsol(ixt,:)=qsol0*Rdefault(index_iso(ixt))
1196      ENDDO
1197#endif       
1198#endif
1199    ENDIF
1200!****************************************************************************************
1201
1202!****************************************************************************************
1203! 2) Initialization to zero
1204!****************************************************************************************
1205!
1206! 2a) Initialization of all argument variables with INTENT(OUT)
1207!****************************************************************************************
1208 cdragh(:)=0. ; cdragm(:)=0.
1209 zu1(:)=0. ; zv1(:)=0.
1210 yus0(:)=0. ; yvs0(:)=0.
1211!albedo SB >>>
1212  alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0.
1213!albedo SB <<<
1214 zxsens(:)=0. ; zxevap(:)=0. ; zxtsol(:)=0. ; zxsnowerosion(:)=0.
1215 d_t_w(:,:)=0. ; d_q_w(:,:)=0. ; d_t_x(:,:)=0. ; d_q_x(:,:)=0.
1216 zxfluxlat(:)=0.
1217 zt2m(:)=0. ; zq2m(:)=0. ; qsat2m(:)=0. ; rh2m(:)=0.
1218 zn2mout(:,:)=0 ;
1219 d_t(:,:)=0. ; d_t_diss(:,:)=0. ; d_q(:,:)=0. ; d_qbs(:,:)=0. ; d_u(:,:)=0. ; d_v(:,:)=0.
1220 zcoefh(:,:,:)=0. ; zcoefm(:,:,:)=0.
1221 zxsens_x(:)=0. ; zxsens_w(:)=0. ; zxfluxlat_x(:)=0. ; zxfluxlat_w(:)=0.
1222 cdragh_x(:)=0. ; cdragh_w(:)=0. ; cdragm_x(:)=0. ; cdragm_w(:)=0.
1223 kh(:)=0. ; kh_x(:)=0. ; kh_w(:)=0.
1224 slab_wfbils(:)=0.
1225 s_pblh(:)=0. ; s_pblh_x(:)=0. ; s_pblh_w(:)=0.
1226 s_plcl(:)=0. ; s_plcl_x(:)=0. ; s_plcl_w(:)=0.
1227 s_capCL(:)=0. ; s_oliqCL(:)=0. ; s_cteiCL(:)=0. ; s_pblT(:)=0.
1228 s_therm(:)=0.
1229 s_trmb1(:)=0. ; s_trmb2(:)=0. ; s_trmb3(:)=0.
1230 zustar(:)=0.
1231 zu10m(:)=0. ; zv10m(:)=0.
1232 fder_print(:)=0.
1233 zxqsurf(:)=0.
1234 delta_qsurf(:) = 0.
1235 zxfluxu(:,:)=0. ; zxfluxv(:,:)=0.
1236 solsw(:,:)=0. ; sollw(:,:)=0.
1237 d_ts(:,:)=0.
1238 evap(:,:)=0.
1239 snowerosion(:,:)=0.
1240 fluxlat(:,:)=0.
1241 wfbils(:,:)=0. ; wfevap(:,:)=0. ;
1242 flux_t(:,:,:)=0. ; flux_q(:,:,:)=0. ; flux_u(:,:,:)=0. ; flux_v(:,:,:)=0.
1243 flux_qbs(:,:,:)=0.
1244 dflux_t(:)=0. ; dflux_q(:)=0.
1245 zxsnow(:)=0.
1246 zxfluxt(:,:)=0. ; zxfluxq(:,:)=0.; zxfluxqbs(:,:)=0.
1247 qsnow(:)=0. ; snowhgt(:)=0. ; to_ice(:)=0. ; sissnow(:)=0.
1248 runoff(:)=0.
1249#ifdef ISO
1250zxxtevap(:,:)=0.
1251 d_xt(:,:,:)=0.
1252 d_xt_x(:,:,:)=0.
1253 d_xt_w(:,:,:)=0.
1254 flux_xt(:,:,:,:)=0.
1255! xtsnow(:,:,:)=0.! attention, xtsnow est l'équivalent de snow et non de qsnow
1256 xtevap(:,:,:)=0.
1257#endif
1258    IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
1259       zcoefh(:,:,:) = 0.0
1260       zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
1261       zcoefm(:,:,:) = 0.0
1262       zcoefm(:,1,:) = 999999. !
1263    ELSE
1264      zcoefm(:,:,is_ave)=0.
1265      zcoefh(:,:,is_ave)=0.
1266    ENDIF
1267!!
1268!  The components "is_ave" of tke_x and wake_deltke are "OUT" variables
1269!jyg<
1270!!    tke(:,:,is_ave)=0.
1271    tke_x(:,:,is_ave)=0.
1272    eps_x(:,:,is_ave)=0.
1273
1274    wake_dltke(:,:,is_ave)=0.
1275!>jyg
1276!!! jyg le 23/02/2013
1277    t2m(:,:)       = 999999.     ! t2m and q2m are meaningfull only over sub-surfaces
1278    q2m(:,:)       = 999999.     ! actually present in the grid cell.
1279!!!
1280    rh2m(:) = 0. ; qsat2m(:) = 0.
1281!!!
1282!!! jyg le 10/02/2012
1283    rh2m_x(:) = 0. ; qsat2m_x(:) = 0. ; rh2m_w(:) = 0. ; qsat2m_w(:) = 0.
1284
1285! 2b) Initialization of all local variables that will be compressed later
1286!****************************************************************************************
1287!!    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
1288    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0
1289!!    zv1 = 0.0     ; yqsurf = 0.0
1290!albedo SB >>>
1291    yqsurf = 0.0  ; yalb = 0.0 ; yalb_vis = 0.0
1292!albedo SB <<<
1293    yrain_f = 0.0 ; ysnow_f = 0.0  ; ybs_f=0.0  ; yfder = 0.0     ; ysolsw = 0.0
1294    ysollw = 0.0  ; yz0m = 0.0 ; yz0h = 0.0    ; yz0h_oupas = 0.0 ; yu1 = 0.0   
1295    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0     ; yqbs1 = 0.0
1296    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
1297    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0
1298    yqbs(:,:)=0.0 
1299    yrugoro = 0.0 ; ywindsp = 0.0   
1300!!    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
1301    yfluxlat=0.0 ; y_flux0(:)=0.0
1302!!    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
1303!!    d_t_diss= 0.0 ;d_u = 0.0     ; d_v = 0.0
1304    yqsol = 0.0   
1305
1306    ytke=0.
1307    yeps=0.
1308    yri0(:)=0.
1309!FC
1310    y_treedrg=0.
1311
1312    ! Martin
1313    ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
1314    yalb3_new = 0.0  ; ysissnow = 0.0
1315    ycldt = 0.0      ; yrmu0 = 0.0
1316    ! Martin
1317    y_d_qbs(:,:)=0.0
1318
1319!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
1320    ytke_x=0.     ; ytke_w=0.        ; ywake_dltke=0.
1321    yeps_x=0.     ; yeps_w=0.
1322    y_d_t_x=0.    ; y_d_t_w=0.       ; y_d_q_x=0.      ; y_d_q_w=0.
1323!!    d_t_w=0.      ; d_q_w=0.         
1324!!    d_t_x=0.      ; d_q_x=0.
1325!!    d_wake_dlt=0.    ; d_wake_dlq=0.
1326    yfluxlat_x=0. ; yfluxlat_w=0.
1327    ywake_s=0.    ; ywake_cstar=0.   ;ywake_dens=0.
1328!!!
1329!!! nrlmd le 13/06/2011
1330    tau_eq=0.     ; delta_coef=0.
1331    y_delta_flux_t1=0.
1332    ydtsurf_th=0.
1333    yts_x(:)=0.      ; yts_w(:)=0.
1334    y_delta_tsurf(:)=0. ; y_delta_qsurf(:)=0.
1335    yqsurf_x(:)=0.      ; yqsurf_w(:)=0.
1336    yg_T(:) = 0. ;        yg_Q(:) = 0.
1337    yGamma_dTs_phiT(:) = 0. ; yGamma_dQs_phiQ(:) = 0.
1338    ydTs_ins(:) = 0. ; ydqs_ins(:) = 0.
1339
1340!!!
1341    ytsoil = 999999.
1342!FC
1343    y_d_u_frein(:,:)=0.
1344    y_d_v_frein(:,:)=0.
1345!FC
1346
1347#ifdef ISO
1348   yxtrain_f = 0.0 ; yxtsnow_f = 0.0
1349   yxtsnow  = 0.0
1350   yxt = 0.0
1351   yxtsol = 0.0
1352   flux_xt = 0.0
1353   yRland_ice = 0.0
1354!   d_xt = 0.0     
1355   y_dflux_xt = 0.0 
1356   dflux_xt=0.0
1357   y_d_xt_x=0.      ; y_d_xt_w=0.       
1358#endif
1359
1360! >> PC
1361!the yfields_out variable is defined in (klon,nbcf_out) even if it is used on
1362!the ORCHIDEE grid and as such should be defined in yfields_out(knon,nbcf_out) but
1363!the knon variable is not known at that level of pbl_surface_mod
1364
1365!the yfields_in variable is defined in (klon,nbcf_in) even if it is used on the
1366!ORCHIDEE grid and as such should be defined in yfields_in(knon,nbcf_in) but the
1367!knon variable is not known at that level of pbl_surface_mod
1368
1369   yfields_out(:,:) = 0.
1370! << PC
1371
1372
1373! 2c) Initialization of all local variables computed within the subsurface loop and used later on
1374!****************************************************************************************
1375    d_t_diss_x(:,:) = 0. ;        d_t_diss_w(:,:) = 0.
1376    d_u_x(:,:)=0. ;               d_u_w(:,:)=0.
1377    d_v_x(:,:)=0. ;               d_v_w(:,:)=0.
1378    flux_t_x(:,:,:)=0. ;          flux_t_w(:,:,:)=0.
1379    flux_q_x(:,:,:)=0. ;          flux_q_w(:,:,:)=0.
1380!
1381!jyg<
1382    flux_u_x(:,:,:)=0. ;          flux_u_w(:,:,:)=0.
1383    flux_v_x(:,:,:)=0. ;          flux_v_w(:,:,:)=0.
1384    fluxlat_x(:,:)=0. ;           fluxlat_w(:,:)=0.
1385!>jyg
1386#ifdef ISO
1387    flux_xt_x(:,:,:,:)=0. ;          flux_xt_w(:,:,:,:)=0.
1388#endif
1389!
1390!jyg<
1391! pblh,plcl,capCL,cteiCL ... are meaningfull only over sub-surfaces
1392! actually present in the grid cell  ==> value set to 999999.
1393!                           
1394!jyg<
1395       ustar(:,:)   = 999999.
1396       wstar(:,:)   = 999999.
1397       windsp(:,:)  = SQRT(u10m(:,:)**2 + v10m(:,:)**2 )
1398       u10m(:,:)    = 999999.
1399       v10m(:,:)    = 999999.
1400!>jyg
1401!
1402       pblh(:,:)   = 999999.        ! Hauteur de couche limite
1403       plcl(:,:)   = 999999.        ! Niveau de condensation de la CLA
1404       capCL(:,:)  = 999999.        ! CAPE de couche limite
1405       oliqCL(:,:) = 999999.        ! eau_liqu integree de couche limite
1406       cteiCL(:,:) = 999999.        ! cloud top instab. crit. couche limite
1407       pblt(:,:)   = 999999.        ! T a la Hauteur de couche limite
1408       therm(:,:)  = 999999.
1409       trmb1(:,:)  = 999999.        ! deep_cape
1410       trmb2(:,:)  = 999999.        ! inhibition
1411       trmb3(:,:)  = 999999.        ! Point Omega
1412!
1413       t2m_x(:,:)    = 999999.
1414       q2m_x(:,:)    = 999999.
1415       ustar_x(:,:)   = 999999.
1416       wstar_x(:,:)   = 999999.
1417       u10m_x(:,:)   = 999999.
1418       v10m_x(:,:)   = 999999.
1419!                           
1420       pblh_x(:,:)   = 999999.      ! Hauteur de couche limite
1421       plcl_x(:,:)   = 999999.      ! Niveau de condensation de la CLA
1422       capCL_x(:,:)  = 999999.      ! CAPE de couche limite
1423       oliqCL_x(:,:) = 999999.      ! eau_liqu integree de couche limite
1424       cteiCL_x(:,:) = 999999.      ! cloud top instab. crit. couche limite
1425       pblt_x(:,:)   = 999999.      ! T a la Hauteur de couche limite
1426       therm_x(:,:)  = 999999.     
1427       trmb1_x(:,:)  = 999999.      ! deep_cape
1428       trmb2_x(:,:)  = 999999.      ! inhibition
1429       trmb3_x(:,:)  = 999999.      ! Point Omega
1430!
1431       t2m_w(:,:)    = 999999.
1432       q2m_w(:,:)    = 999999.
1433       ustar_w(:,:)   = 999999.
1434       wstar_w(:,:)   = 999999.
1435       u10m_w(:,:)   = 999999.
1436       v10m_w(:,:)   = 999999.
1437                           
1438       pblh_w(:,:)   = 999999.      ! Hauteur de couche limite
1439       plcl_w(:,:)   = 999999.      ! Niveau de condensation de la CLA
1440       capCL_w(:,:)  = 999999.      ! CAPE de couche limite
1441       oliqCL_w(:,:) = 999999.      ! eau_liqu integree de couche limite
1442       cteiCL_w(:,:) = 999999.      ! cloud top instab. crit. couche limite
1443       pblt_w(:,:)   = 999999.      ! T a la Hauteur de couche limite
1444       therm_w(:,:)  = 999999.     
1445       trmb1_w(:,:)  = 999999.      ! deep_cape
1446       trmb2_w(:,:)  = 999999.      ! inhibition
1447       trmb3_w(:,:)  = 999999.      ! Point Omega
1448!!!     
1449!
1450!!!
1451!****************************************************************************************
1452! 3) - Calculate pressure thickness of each layer
1453!    - Calculate the wind at first layer
1454!    - Mean calculations of albedo
1455!    - Calculate net radiance at sub-surface
1456!****************************************************************************************
1457    DO k = 1, klev
1458       DO i = 1, klon
1459          delp(i,k) = paprs(i,k)-paprs(i,k+1)
1460       ENDDO
1461    ENDDO
1462
1463!****************************************************************************************
1464! Test for rugos........ from physiq.. A la fin plutot???
1465!
1466!****************************************************************************************
1467
1468    DO nsrf = 1, nbsrf
1469       DO i = 1, klon
1470          z0m(i,nsrf) = MAX(z0m(i,nsrf),z0min)
1471          z0h(i,nsrf) = MAX(z0h(i,nsrf),z0min)
1472       ENDDO
1473    ENDDO
1474
1475! Mean calculations of albedo
1476!
1477! * alb  : mean albedo for whole SW interval
1478!
1479! Mean albedo for grid point
1480! * alb_m  : mean albedo at whole SW interval
1481
1482    alb_dir_m(:,:) = 0.0
1483    alb_dif_m(:,:) = 0.0
1484    DO k = 1, nsw
1485     DO nsrf = 1, nbsrf
1486       DO i = 1, klon
1487          alb_dir_m(i,k) = alb_dir_m(i,k) + alb_dir(i,k,nsrf) * pctsrf(i,nsrf)
1488          alb_dif_m(i,k) = alb_dif_m(i,k) + alb_dif(i,k,nsrf) * pctsrf(i,nsrf)
1489       ENDDO
1490     ENDDO
1491    ENDDO
1492
1493! We here suppose the fraction f1 of incoming radiance of visible radiance
1494! as a fraction of all shortwave radiance
1495    f1 = 0.5
1496!    f1 = 1    ! put f1=1 to recreate old calculations
1497
1498!f1 is already included with SFRWL values in each surf files
1499    alb=0.0
1500    DO k=1,nsw
1501      DO nsrf = 1, nbsrf
1502        DO i = 1, klon
1503            alb(i,nsrf) = alb(i,nsrf) + alb_dir(i,k,nsrf)*SFRWL(k)
1504        ENDDO
1505      ENDDO
1506    ENDDO
1507
1508    alb_m=0.0
1509    DO k = 1,nsw
1510      DO i = 1, klon
1511        alb_m(i) = alb_m(i) + alb_dir_m(i,k)*SFRWL(k)
1512      ENDDO
1513    ENDDO
1514!albedo SB <<<
1515
1516
1517
1518! Calculation of mean temperature at surface grid points
1519    ztsol(:) = 0.0
1520    DO nsrf = 1, nbsrf
1521       DO i = 1, klon
1522          ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
1523       ENDDO
1524    ENDDO
1525
1526! Linear distrubution on sub-surface of long- and shortwave net radiance
1527    DO nsrf = 1, nbsrf
1528       DO i = 1, klon
1529          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
1530!--OB this line is not satisfactory because alb is the direct albedo not total albedo
1531          solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
1532       ENDDO
1533    ENDDO
1534!
1535!<al1: second order corrections
1536!- net = dwn -up; up=sig( T4 + 4sum%T3T' + 6sum%T2T'2 +...)
1537   IF (iflag_order2_sollw == 1) THEN
1538    meansqT(:) = 0. ! as working buffer
1539    DO nsrf = 1, nbsrf
1540     DO i = 1, klon
1541      meansqT(i) = meansqT(i)+(ts(i,nsrf)-ztsol(i))**2 *pctsrf(i,nsrf)
1542     ENDDO
1543    ENDDO
1544    DO nsrf = 1, nbsrf
1545     DO i = 1, klon
1546      sollw(i,nsrf) = sollw(i,nsrf) &
1547                + 6.0*RSIGMA*ztsol(i)**2 *(meansqT(i)-(ztsol(i)-ts(i,nsrf))**2)
1548     ENDDO
1549    ENDDO
1550   ENDIF   ! iflag_order2_sollw == 1
1551!>al1
1552
1553!--OB add diffuse fraction of SW down
1554   DO n=1,nbcf_out
1555       IF (cfname_out(n) == "swdownfdiff" ) fields_out(:,n) = solswfdiff_m(:)
1556   ENDDO
1557! >> PC
1558   IF (carbon_cycle_cpl .AND. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
1559       r_co2_ppm(:) = co2_send(:)
1560       DO n=1,nbcf_out
1561           IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_send(:)
1562       ENDDO
1563   ENDIF
1564   IF ( .NOT. carbon_cycle_tr .AND. nbcf_out.GT.0 ) THEN
1565       r_co2_ppm(:) = co2_ppm     ! Constant field
1566       DO n=1,nbcf_out
1567           IF (cfname_out(n) == "atmco2" ) fields_out(:,n) = co2_ppm
1568       ENDDO
1569   ENDIF
1570! << PC
1571
1572!****************************************************************************************
1573! 4) Loop over different surfaces
1574!
1575! Only points containing a fraction of the sub surface will be treated.
1576!
1577!****************************************************************************************
1578                                                                          !<<<<<<<<<<<<<
1579    loop_nbsrf: DO nsrf = 1, nbsrf                                        !<<<<<<<<<<<<<
1580                                                                          !<<<<<<<<<<<<<
1581       IF (prt_level >=10) print *,' Loop nsrf ',nsrf
1582!
1583       IF (iflag_split_ref == 3) THEN
1584         IF (nsrf == is_oce) THEN
1585            iflag_split = 1
1586         ELSE
1587            iflag_split=0
1588         ENDIF   !! (nsrf == is_oce)
1589       ELSE                     
1590         iflag_split = iflag_split_ref
1591       ENDIF   !! (iflag_split_ref == 3)
1592
1593! Search for index(ni) and size(knon) of domaine to treat
1594       ni(:) = 0
1595       knon  = 0
1596       DO i = 1, klon
1597          IF (pctsrf(i,nsrf) > 0.) THEN
1598             knon = knon + 1
1599             ni(knon) = i
1600          ENDIF
1601       ENDDO
1602
1603!!! jyg le 19/08/2012
1604!       IF (knon <= 0) THEN
1605!         IF (prt_level >= 10) print *,' no grid point for nsrf= ',nsrf
1606!         cycle loop_nbsrf
1607!       ENDIF
1608!!!
1609
1610       ! write index, with IOIPSL
1611       IF (debugindex .AND. mpi_size==1) THEN
1612          tabindx(:)=0.
1613          DO i=1,knon
1614             tabindx(i)=REAL(i)
1615          ENDDO
1616          debugtab(:,:) = 0.
1617          ndexbg(:) = 0
1618          CALL gath2cpl(tabindx,debugtab,knon,ni)
1619          CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,nbp_lon*nbp_lat, ndexbg)
1620       ENDIF
1621       
1622!****************************************************************************************
1623! 5) Compress variables
1624!
1625!****************************************************************************************
1626
1627!
1628!jyg<    (20190926)
1629!   Provisional : set ybeta to standard values
1630       IF (nsrf .NE. is_ter) THEN
1631           ybeta(1:knon) = 1.
1632       ELSE
1633           IF (iflag_split .EQ. 0) THEN
1634              ybeta(1:knon) = 1.
1635           ELSE
1636             DO j = 1, knon
1637                i = ni(j)
1638                ybeta(j)   = beta(i,nsrf)
1639             ENDDO
1640           ENDIF  ! (iflag_split .LE.1)
1641       ENDIF !  (nsrf .NE. is_ter)
1642!>jyg
1643!
1644       DO j = 1, knon
1645          i = ni(j)
1646          ypct(j)    = pctsrf(i,nsrf)
1647          yts(j)     = ts(i,nsrf)
1648          ysnow(j)   = snow(i,nsrf)
1649          yqsurf(j)  = qsurf(i,nsrf)
1650          yalb(j)    = alb(i,nsrf)
1651!albedo SB >>>
1652          yalb_vis(j) = alb_dir(i,1,nsrf)
1653          IF (nsw==6) THEN
1654            yalb_vis(j)=(alb_dir(i,1,nsrf)*SFRWL(1)+alb_dir(i,2,nsrf)*SFRWL(2) &
1655              +alb_dir(i,3,nsrf)*SFRWL(3))/(SFRWL(1)+SFRWL(2)+SFRWL(3))
1656          ENDIF
1657!albedo SB <<<
1658          yrain_f(j) = rain_f(i)
1659          ysnow_f(j) = snow_f(i)
1660          ybs_f(j)   = bs_f(i)
1661          yagesno(j) = agesno(i,nsrf)
1662          yfder(j)   = fder(i)
1663          ylwdown(j) = lwdown_m(i)
1664          ygustiness(j) = gustiness(i)
1665          ysolsw(j)  = solsw(i,nsrf)
1666          ysollw(j)  = sollw(i,nsrf)
1667          yz0m(j)  = z0m(i,nsrf)
1668          yz0h(j)  = z0h(i,nsrf)
1669          yrugoro(j) = rugoro(i)
1670          yu1(j)     = u(i,1)
1671          yv1(j)     = v(i,1)
1672          yqbs1(j)   = qbs(i,1)
1673          ypaprs(j,klev+1) = paprs(i,klev+1)
1674!jyg<
1675!!          ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )
1676          ywindsp(j) = windsp(i,nsrf)
1677!>jyg
1678          ! Martin and Etienne
1679          yzmea(j)   = zmea(i)
1680          yzsig(j)   = zsig(i)
1681          ycldt(j)   = cldt(i)
1682          yrmu0(j)   = rmu0(i)
1683          ! Martin
1684!!! nrlmd le 13/06/2011
1685          y_delta_tsurf(j)=delta_tsurf(i,nsrf)
1686          yfluxbs(j)=0.0
1687          y_flux_bs(j) = 0.0
1688!!!
1689#ifdef ISO
1690          DO ixt=1,ntraciso
1691            yxtrain_f(ixt,j) = xtrain_f(ixt,i)
1692            yxtsnow_f(ixt,j) = xtsnow_f(ixt,i) 
1693          ENDDO
1694          DO ixt=1,niso
1695            yxtsnow(ixt,j)   = xtsnow(ixt,i,nsrf)
1696          ENDDO   
1697          !IF (nsrf == is_lic) THEN
1698          DO ixt=1,niso
1699            yRland_ice(ixt,j)= Rland_ice(ixt,i) 
1700          ENDDO   
1701          !endif !IF (nsrf == is_lic) THEN
1702#ifdef ISOVERIF
1703          IF (iso_eau >= 0) THEN
1704              call iso_verif_egalite_choix(ysnow_f(j), &
1705     &          yxtsnow_f(iso_eau,j),'pbl_surf_mod 862', &
1706     &          errmax,errmaxrel)
1707              call iso_verif_egalite_choix(ysnow(j), &
1708     &          yxtsnow(iso_eau,j),'pbl_surf_mod 872', &
1709     &          errmax,errmaxrel)
1710          ENDIF
1711#endif
1712#ifdef ISOVERIF
1713         DO ixt=1,ntraciso
1714           call iso_verif_noNaN(yxtsnow_f(ixt,j),'pbl_surf_mod 921')
1715         ENDDO
1716#endif
1717#endif
1718       ENDDO
1719! >> PC
1720!--compressing fields_out onto ORCHIDEE grid
1721!--these fields are shared and used directly surf_land_orchidee_mod
1722       DO n = 1, nbcf_out
1723         DO j = 1, knon
1724           i = ni(j)
1725           yfields_out(j,n) = fields_out(i,n)
1726         ENDDO
1727       ENDDO
1728! << PC
1729       DO k = 1, klev
1730          DO j = 1, knon
1731             i = ni(j)
1732             ypaprs(j,k) = paprs(i,k)
1733             ypplay(j,k) = pplay(i,k)
1734             ydelp(j,k)  = delp(i,k)
1735          ENDDO
1736       ENDDO
1737!
1738!!! jyg le 07/02/2012 et le 10/04/2013
1739        DO k = 1, klev+1
1740          DO j = 1, knon
1741             i = ni(j)
1742!jyg<
1743!!             ytke(j,k)   = tke(i,k,nsrf)
1744             ytke(j,k)   = tke_x(i,k,nsrf)
1745          ENDDO
1746        ENDDO
1747!>jyg
1748        DO k = 1, klev
1749          DO j = 1, knon
1750             i = ni(j)
1751             y_treedrg(j,k) =  treedrg(i,k,nsrf)
1752             yu(j,k) = u(i,k)
1753             yv(j,k) = v(i,k)
1754             yt(j,k) = t(i,k)
1755             yq(j,k) = q(i,k)
1756             yqbs(j,k)=qbs(i,k)
1757#ifdef ISO
1758             DO ixt=1,ntraciso   
1759               yxt(ixt,j,k) = xt(ixt,i,k)
1760             ENDDO !DO ixt=1,ntraciso
1761#endif
1762          ENDDO
1763        ENDDO
1764!
1765       IF (iflag_split.GE.1) THEN
1766!!! nrlmd le 02/05/2011
1767        DO k = 1, klev
1768          DO j = 1, knon
1769             i = ni(j)
1770             yu_x(j,k) = u(i,k)
1771             yv_x(j,k) = v(i,k)
1772             yt_x(j,k) = t(i,k)-wake_s(i)*wake_dlt(i,k)
1773             yq_x(j,k) = q(i,k)-wake_s(i)*wake_dlq(i,k)
1774             yu_w(j,k) = u(i,k)
1775             yv_w(j,k) = v(i,k)
1776             yt_w(j,k) = t(i,k)+(1.-wake_s(i))*wake_dlt(i,k)
1777             yq_w(j,k) = q(i,k)+(1.-wake_s(i))*wake_dlq(i,k)
1778!!!
1779#ifdef ISO
1780             DO ixt=1,ntraciso
1781               yxt_x(ixt,j,k) = xt(ixt,i,k)-wake_s(i)*wake_dlxt(ixt,i,k)
1782               yxt_w(ixt,j,k) = xt(ixt,i,k)+(1.-wake_s(i))*wake_dlxt(ixt,i,k)
1783             ENDDO
1784#endif
1785          ENDDO
1786        ENDDO
1787
1788        IF (prt_level .ge. 10) THEN
1789          print *,'pbl_surface, wake_s(1), wake_dlt(1,:) ', wake_s(1), wake_dlt(1,:)
1790          print *,'pbl_surface, wake_s(1), wake_dlq(1,:) ', wake_s(1), wake_dlq(1,:)
1791        ENDIF
1792
1793!!! nrlmd le 02/05/2011
1794        DO k = 1, klev+1
1795          DO j = 1, knon
1796             i = ni(j)
1797!jyg<
1798!!             ytke_x(j,k) = tke(i,k,nsrf)-wake_s(i)*wake_dltke(i,k,nsrf)
1799!!             ytke_w(j,k) = tke(i,k,nsrf)+(1.-wake_s(i))*wake_dltke(i,k,nsrf)
1800!!             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
1801!!             ytke(j,k)     = tke(i,k,nsrf)
1802!
1803             ytke_x(j,k)      = tke_x(i,k,nsrf)
1804             ytke(j,k)        = tke_x(i,k,nsrf)+wake_s(i)*wake_dltke(i,k,nsrf)
1805             ytke_w(j,k)      = tke_x(i,k,nsrf)+wake_dltke(i,k,nsrf)
1806             ywake_dltke(j,k) = wake_dltke(i,k,nsrf)
1807           
1808!>jyg
1809          ENDDO
1810        ENDDO
1811!!!
1812!!! jyg le 07/02/2012
1813        DO j = 1, knon
1814          i = ni(j)
1815          ywake_s(j)=wake_s(i)
1816          ywake_cstar(j)=wake_cstar(i)
1817          ywake_dens(j)=wake_dens(i)
1818        ENDDO
1819!!!
1820!!! nrlmd le 13/06/2011
1821        DO j=1,knon
1822         yts_x(j)=yts(j)-ywake_s(j)*y_delta_tsurf(j)
1823         yts_w(j)=yts(j)+(1.-ywake_s(j))*y_delta_tsurf(j)
1824        ENDDO
1825!!!
1826       ENDIF  ! (iflag_split .ge.1)
1827!!!
1828       DO k = 1, nsoilmx
1829          DO j = 1, knon
1830             i = ni(j)
1831             ytsoil(j,k) = ftsoil(i,k,nsrf)
1832          ENDDO
1833       ENDDO
1834       
1835       ! qsol(water height in soil) only for bucket continental model
1836       IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
1837          DO j = 1, knon
1838             i = ni(j)
1839             yqsol(j) = qsol(i)
1840#ifdef ISO
1841             DO ixt=1,niso
1842               yxtsol(ixt,j) = xtsol(ixt,i)
1843             ENDDO
1844#endif
1845          ENDDO
1846       ENDIF
1847
1848       if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
1849          if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
1850             ydelta_sal(:knon) = delta_sal(ni(:knon))
1851             ydelta_sst(:knon) = delta_sst(ni(:knon))
1852             ydter(:knon) = dter(ni(:knon))
1853             ydser(:knon) = dser(ni(:knon))
1854             ydt_ds(:knon) = dt_ds(ni(:knon))
1855          end if
1856         
1857          yds_ns(:knon) = ds_ns(ni(:knon))
1858          ydt_ns(:knon) = dt_ns(ni(:knon))
1859       end if
1860       
1861!****************************************************************************************
1862! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
1863!
1864!****************************************************************************************
1865
1866
1867!!! jyg le 07/02/2012
1868       IF (iflag_split .eq.0) THEN
1869!!!
1870!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1871! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1872!       CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1873!           yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
1874!           yts, yqsurf, yrugos, &
1875!           ycdragm, ycdragh )
1876! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1877        DO i = 1, knon
1878!          print*,'PBL ',i,RD
1879!          print*,'PBL ',yt(i,1),ypaprs(i,1),ypplay(i,1)
1880           zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1881                * (ypaprs(i,1)-ypplay(i,1))
1882           speed(i) = SQRT(yu(i,1)**2+yv(i,1)**2)
1883        ENDDO
1884        CALL cdrag(knon, nsrf, &
1885            speed, yt(:,1), yq(:,1), zgeo1, ypaprs(:,1), s_pblh, &
1886            yts, yqsurf, yz0m, yz0h, yri0, 0, &
1887            ycdragm, ycdragh, zri1, pref, rain_f, zxtsol, ypplay(:,1))
1888
1889! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
1890     IF (ok_prescr_ust) THEN
1891      DO i = 1, knon
1892       print *,'ycdragm avant=',ycdragm(i)
1893       vent= sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))
1894!      ycdragm(i) = ust*ust/(1.+(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)))
1895!      ycdragm(i) = ust*ust/((1.+sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1))) &
1896!     *sqrt(yu(i,1)*yu(i,1)+yv(i,1)*yv(i,1)))
1897       ycdragm(i) = ust*ust/(1.+vent)/vent
1898!      print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1)
1899      ENDDO
1900     ENDIF
1901
1902        IF (prt_level >=10) print *,'cdrag -> ycdragh ', ycdragh(1:knon)
1903       ELSE  !(iflag_split .eq.0)
1904
1905! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1906!       CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1907!           yu_x(:,1), yv_x(:,1), yt_x(:,1), yq_x(:,1), &
1908!           yts_x, yqsurf, yrugos, &
1909!           ycdragm_x, ycdragh_x )
1910! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1911        DO i = 1, knon
1912           zgeo1_x(i) = RD * yt_x(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1913                * (ypaprs(i,1)-ypplay(i,1))
1914           speed_x(i) = SQRT(yu_x(i,1)**2+yv_x(i,1)**2)
1915        ENDDO
1916
1917
1918            CALL cdrag(knon, nsrf, &
1919            speed_x, yt_x(:,1), yq_x(:,1), zgeo1_x, ypaprs(:,1),s_pblh_x,&
1920            yts_x, yqsurf_x, yz0m, yz0h, yri0, 0, &
1921            ycdragm_x, ycdragh_x, zri1_x, pref_x, rain_f, zxtsol, ypplay(:,1) )
1922
1923! --- special Dice. JYG+MPL 25112013
1924        IF (ok_prescr_ust) THEN
1925         DO i = 1, knon
1926!         print *,'ycdragm_x avant=',ycdragm_x(i)
1927          vent= sqrt(yu_x(i,1)*yu_x(i,1)+yv_x(i,1)*yv_x(i,1))
1928          ycdragm_x(i) = ust*ust/(1.+vent)/vent
1929!         print *,'ycdragm_x ust yu yv apres=',ycdragm_x(i),ust,yu_x(i,1),yv_x(i,1)
1930         ENDDO
1931        ENDIF
1932        IF (prt_level >=10) print *,'clcdrag -> ycdragh_x ', ycdragh_x(1:knon)
1933!
1934! Faire disparaitre les lignes commentees fin 2015 (le temps des tests)
1935!        CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
1936!            yu_w(:,1), yv_w(:,1), yt_w(:,1), yq_w(:,1), &
1937!            yts_w, yqsurf, yz0m, &
1938!            ycdragm_w, ycdragh_w )
1939! Fuxing WANG, 04/03/2015, replace the clcdrag by the merged version: cdrag
1940        DO i = 1, knon
1941           zgeo1_w(i) = RD * yt_w(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
1942                * (ypaprs(i,1)-ypplay(i,1))
1943           speed_w(i) = SQRT(yu_w(i,1)**2+yv_w(i,1)**2)
1944        ENDDO
1945        CALL cdrag(knon, nsrf, &
1946            speed_w, yt_w(:,1), yq_w(:,1), zgeo1_w, ypaprs(:,1),s_pblh_w,&
1947            yts_w, yqsurf_w, yz0m, yz0h, yri0, 0, &
1948            ycdragm_w, ycdragh_w, zri1_w, pref_w, rain_f, zxtsol, ypplay(:,1) )
1949!
1950        IF(ok_bug_zg_wk_pbl) THEN
1951         zgeo1(1:knon) = wake_s(1:knon)*zgeo1_w(1:knon) + (1.-wake_s(1:knon))*zgeo1_x(1:knon)
1952        ELSE
1953         zgeo1(1:knon) = ywake_s(1:knon)*zgeo1_w(1:knon) + (1.-ywake_s(1:knon))*zgeo1_x(1:knon)
1954        ENDIF
1955
1956! --- special Dice. JYG+MPL 25112013 puis BOMEX
1957        IF (ok_prescr_ust) THEN
1958         DO i = 1, knon
1959!         print *,'ycdragm_w avant=',ycdragm_w(i)
1960          vent= sqrt(yu_w(i,1)*yu_w(i,1)+yv_w(i,1)*yv_w(i,1))
1961          ycdragm_w(i) = ust*ust/(1.+vent)/vent
1962!         print *,'ycdragm_w ust yu yv apres=',ycdragm_w(i),ust,yu_w(i,1),yv_w(i,1)
1963         ENDDO
1964        ENDIF
1965        IF (prt_level >=10) print *,'clcdrag -> ycdragh_w ', ycdragh_w(1:knon)
1966!!!
1967       ENDIF  ! (iflag_split .eq.0)
1968!!!
1969       
1970
1971!****************************************************************************************
1972! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefh et ycoefm.
1973!
1974!****************************************************************************************
1975
1976!!! jyg le 07/02/2012
1977       IF (iflag_split .eq.0) THEN
1978!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
1979      IF (prt_level >=10) THEN
1980      print *,' args coef_diff_turb: yu ',  yu(1:knon,:) 
1981      print *,' args coef_diff_turb: yv ',  yv(1:knon,:)   
1982      print *,' args coef_diff_turb: yq ',  yq(1:knon,:)   
1983      print *,' args coef_diff_turb: yt ',  yt(1:knon,:)   
1984      print *,' args coef_diff_turb: yts ', yts(1:knon)
1985      print *,' args coef_diff_turb: yz0m ', yz0m(1:knon)
1986      print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon) 
1987      print *,' args coef_diff_turb: ycdragm ', ycdragm(1:knon)
1988      print *,' args coef_diff_turb: ycdragh ', ycdragh(1:knon)
1989      print *,' args coef_diff_turb: ytke ', ytke(1:knon,:)   
1990       ENDIF
1991
1992        IF (iflag_pbl>=50) THEN
1993        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm(1:knon), ycdragh(1:knon),yus0(1:knon),yvs0(1:knon),yts(1:knon), &
1994                  yu(1:knon,:),yv(1:knon,:),yt(1:knon,:),yq(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),       &
1995                  ytke(1:knon,:),yeps(1:knon,:), ycoefm(1:knon,:), ycoefh(1:knon,:))
1996
1997        ELSE
1998
1999        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
2000            ypaprs, ypplay, yu, yv, yq, yt, yts, yqsurf, ycdragm, &
2001            ycoefm, ycoefh, ytke, yeps, y_treedrg)
2002!            ycoefm, ycoefh, ytke)
2003!FC y_treedrg ajoute
2004       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
2005! In this case, coef_diff_turb is called for the Cd only
2006       DO k = 2, klev
2007          DO j = 1, knon
2008             i = ni(j)
2009             ycoefh(j,k)   = zcoefh(i,k,nsrf)
2010             ycoefm(j,k)   = zcoefm(i,k,nsrf)
2011          ENDDO
2012       ENDDO
2013       ENDIF
2014
2015       ENDIF ! iflag_pbl >= 50
2016
2017        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh ',ycoefh(1:knon,:)
2018
2019
2020       ELSE  !(iflag_split .eq.0)
2021
2022     
2023      IF (prt_level >=10) THEN
2024      print *,' args coef_diff_turb: yu_x ',  yu_x(1:knon,:)     
2025      print *,' args coef_diff_turb: yv_x ',  yv_x(1:knon,:)     
2026      print *,' args coef_diff_turb: yq_x ',  yq_x(1:knon,:)     
2027      print *,' args coef_diff_turb: yt_x ',  yt_x(1:knon,:)     
2028      print *,' args coef_diff_turb: yts_x ', yts_x(1:knon)
2029      print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon) 
2030      print *,' args coef_diff_turb: ycdragm_x ', ycdragm_x(1:knon)
2031      print *,' args coef_diff_turb: ycdragh_x ', ycdragh_x(1:knon)
2032      print *,' args coef_diff_turb: ytke_x ', ytke_x(1:knon,:)   
2033      ENDIF
2034
2035
2036        IF (iflag_pbl>=50) THEN
2037     
2038        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_x(1:knon),ycdragh_x(1:knon),yus0(1:knon),yvs0(1:knon),yts_x(1:knon),    &
2039                       yu_x(1:knon,:),yv_x(1:knon,:),yt_x(1:knon,:),yq_x(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),  &
2040                       ytke_x(1:knon,:),yeps_x(1:knon,:),ycoefm_x(1:knon,:), ycoefh_x(1:knon,:))
2041
2042        ELSE
2043
2044        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
2045            ypaprs, ypplay, yu_x, yv_x, yq_x, yt_x, yts_x, yqsurf_x, ycdragm_x, &
2046            ycoefm_x, ycoefh_x, ytke_x,yeps_x,y_treedrg)
2047!            ycoefm_x, ycoefh_x, ytke_x)
2048!FC doit on le mettre ( on ne l utilise pas si il y a du spliting)
2049       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
2050! In this case, coef_diff_turb is called for the Cd only
2051       DO k = 2, klev
2052          DO j = 1, knon
2053             i = ni(j)
2054             ycoefh_x(j,k)   = zcoefh(i,k,nsrf)
2055             ycoefm_x(j,k)   = zcoefm(i,k,nsrf)
2056          ENDDO
2057       ENDDO
2058       ENDIF
2059
2060        ENDIF ! iflag_pbl >= 50
2061
2062        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_x ',ycoefh_x(1:knon,:)
2063!
2064      IF (prt_level >=10) THEN
2065      print *,' args coef_diff_turb: yu_w ',  yu_w(1:knon,:)
2066      print *,' args coef_diff_turb: yv_w ',  yv_w(1:knon,:) 
2067      print *,' args coef_diff_turb: yq_w ',  yq_w(1:knon,:) 
2068      print *,' args coef_diff_turb: yt_w ',  yt_w(1:knon,:) 
2069      print *,' args coef_diff_turb: yts_w ', yts_w(1:knon)
2070      print *,' args coef_diff_turb: yqsurf ', yqsurf(1:knon) 
2071      print *,' args coef_diff_turb: ycdragm_w ', ycdragm_w(1:knon)
2072      print *,' args coef_diff_turb: ycdragh_w ', ycdragh_w(1:knon)
2073      print *,' args coef_diff_turb: ytke_w ', ytke_w(1:knon,:)
2074      ENDIF
2075     
2076        IF (iflag_pbl>=50) THEN
2077       
2078        CALL call_atke(dtime,knon,klev,nsrf,ni,ycdragm_w(1:knon),ycdragh_w(1:knon),yus0(1:knon),yvs0(1:knon),yts_w(1:knon), &
2079                yu_w(1:knon,:),yv_w(1:knon,:),yt_w(1:knon,:),yq_w(1:knon,:),ypplay(1:knon,:),ypaprs(1:knon,:),      &
2080                ytke_w(1:knon,:),yeps_w(1:knon,:),ycoefm_w(1:knon,:),ycoefh_w(1:knon,:))
2081
2082        ELSE
2083
2084        CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
2085            ypaprs, ypplay, yu_w, yv_w, yq_w, yt_w, yts_w, yqsurf_w, ycdragm_w, &
2086            ycoefm_w, ycoefh_w, ytke_w,yeps_w,y_treedrg)
2087!            ycoefm_w, ycoefh_w, ytke_w)
2088       IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
2089! In this case, coef_diff_turb is called for the Cd only
2090       DO k = 2, klev
2091          DO j = 1, knon
2092             i = ni(j)
2093             ycoefh_w(j,k)   = zcoefh(i,k,nsrf)
2094             ycoefm_w(j,k)   = zcoefm(i,k,nsrf)
2095          ENDDO
2096       ENDDO
2097       ENDIF
2098
2099       ENDIF ! iflag_pbl >= 50
2100
2101
2102        IF (prt_level >=10) print *,'coef_diff_turb -> ycoefh_w ',ycoefh_w(1:knon,:)
2103
2104!!!jyg le 10/04/2013
2105!!   En attendant de traiter le transport des traceurs dans les poches froides, formule
2106!!   arbitraire pour ycoefh et ycoefm
2107      DO k = 2,klev
2108        DO j = 1,knon
2109         ycoefh(j,k) = ycoefh_x(j,k) + ywake_s(j)*(ycoefh_w(j,k) - ycoefh_x(j,k))
2110         ycoefm(j,k) = ycoefm_x(j,k) + ywake_s(j)*(ycoefm_w(j,k) - ycoefm_x(j,k))
2111        ENDDO
2112      ENDDO
2113
2114
2115       ENDIF  ! (iflag_split .eq.0)
2116
2117       
2118!****************************************************************************************
2119!
2120! 8) "La descente" - "The downhill"
2121
2122!  climb_hq_down and climb_wind_down calculate the coefficients
2123!  Ccoef_X et Dcoef_X for X=[H, Q, U, V].
2124!  Only the coefficients at surface for H and Q are returned.
2125!
2126!****************************************************************************************
2127
2128! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
2129!!! jyg le 07/02/2012
2130       IF (iflag_split .eq.0) THEN
2131!!!
2132!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2133        CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
2134            ydelp, yt, yq, dtime, &
2135!!! jyg le 09/05/2011
2136            CcoefH, CcoefQ, DcoefH, DcoefQ, &
2137            Kcoef_hq, gama_q, gama_h, &
2138!!!
2139            AcoefH, AcoefQ, BcoefH, BcoefQ &
2140#ifdef ISO
2141         &   ,yxt, CcoefXT, DcoefXT, gama_xt, AcoefXT, BcoefXT &
2142#endif               
2143         &   )
2144       ELSE  !(iflag_split .eq.0)
2145        CALL climb_hq_down(knon, ycoefh_x, ypaprs, ypplay, &
2146            ydelp, yt_x, yq_x, dtime, &
2147!!! nrlmd le 02/05/2011
2148            CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
2149            Kcoef_hq_x, gama_q_x, gama_h_x, &
2150!!!
2151            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x &
2152#ifdef ISO
2153         &   ,yxt_x, CcoefXT_x, DcoefXT_x, gama_xt_x, AcoefXT_x, BcoefXT_x &
2154#endif               
2155         &   )
2156!!!
2157       IF (prt_level >=10) THEN
2158         PRINT *,'pbl_surface (climb_hq_down.x->) AcoefH_x ',AcoefH_x
2159         PRINT *,'pbl_surface (climb_hq_down.x->) AcoefQ_x ',AcoefQ_x
2160         PRINT *,'pbl_surface (climb_hq_down.x->) BcoefH_x ',BcoefH_x
2161         PRINT *,'pbl_surface (climb_hq_down.x->) BcoefQ_x ',BcoefQ_x
2162       ENDIF
2163!
2164        CALL climb_hq_down(knon, ycoefh_w, ypaprs, ypplay, &
2165            ydelp, yt_w, yq_w, dtime, &
2166!!! nrlmd le 02/05/2011
2167            CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
2168            Kcoef_hq_w, gama_q_w, gama_h_w, &
2169!!!
2170            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w &
2171#ifdef ISO
2172         &   ,yxt_w, CcoefXT_w, DcoefXT_w, gama_xt_w, AcoefXT_w, BcoefXT_w &
2173#endif               
2174         &   )
2175!!!
2176       IF (prt_level >=10) THEN
2177         PRINT *,'pbl_surface (climb_hq_down.w->) AcoefH_w ',AcoefH_w
2178         PRINT *,'pbl_surface (climb_hq_down.w->) AcoefQ_w ',AcoefQ_w
2179         PRINT *,'pbl_surface (climb_hq_down.w->) BcoefH_w ',BcoefH_w
2180         PRINT *,'pbl_surface (climb_hq_down.w->) BcoefQ_w ',BcoefQ_w
2181       ENDIF
2182!!!
2183       ENDIF  ! (iflag_split .eq.0)
2184!!!
2185
2186! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
2187!!! jyg le 07/02/2012
2188       IF (iflag_split .eq.0) THEN
2189!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2190        CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
2191!!! jyg le 09/05/2011
2192            CcoefU, CcoefV, DcoefU, DcoefV, &
2193            Kcoef_m, alf_1, alf_2, &
2194!!!
2195            AcoefU, AcoefV, BcoefU, BcoefV)
2196       ELSE  ! (iflag_split .eq.0)
2197        CALL climb_wind_down(knon, dtime, ycoefm_x, ypplay, ypaprs, yt_x, ydelp, yu_x, yv_x, &
2198!!! nrlmd le 02/05/2011
2199            CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
2200            Kcoef_m_x, alf_1_x, alf_2_x, &
2201!!!
2202            AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x)
2203!
2204        CALL climb_wind_down(knon, dtime, ycoefm_w, ypplay, ypaprs, yt_w, ydelp, yu_w, yv_w, &
2205!!! nrlmd le 02/05/2011
2206            CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
2207            Kcoef_m_w, alf_1_w, alf_2_w, &
2208!!!
2209            AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w)
2210!!!     
2211       ENDIF  ! (iflag_split .eq.0)
2212!!!
2213
2214! For blowing snow:
2215    IF (ok_bs) THEN
2216     ! following Bintanja et al 2000, part II and Vionnet V PhD thesis
2217     ! we assume that the eddy diffsivity coefficient for
2218     ! suspended particles is a fraction of Kh
2219     do k=1,klev
2220        do j=1,knon
2221           ycoefqbs(j,k)=ycoefh(j,k)*zeta_bs
2222        enddo
2223     enddo
2224     CALL climb_qbs_down(knon, ycoefqbs, ypaprs, ypplay, &
2225     ydelp, yt, yqbs, dtime, &
2226     CcoefQBS, DcoefQBS, &
2227     Kcoef_qbs, gama_qbs, &
2228     AcoefQBS, BcoefQBS)
2229    ENDIF
2230
2231!****************************************************************************************
2232! 9) Small calculations
2233!
2234!****************************************************************************************
2235
2236! - Reference pressure is given the values at surface level         
2237       ypsref(:) = ypaprs(:,1) 
2238
2239! - CO2 field on 2D grid to be sent to ORCHIDEE
2240!   Transform to compressed field
2241       IF (carbon_cycle_cpl) THEN
2242          DO i=1,knon
2243             r_co2_ppm(i) = co2_send(ni(i))
2244          ENDDO
2245       ELSE
2246          r_co2_ppm(:) = co2_ppm     ! Constant field
2247       ENDIF
2248
2249!!! nrlmd le 02/05/2011  -----------------------On raccorde les 2 colonnes dans la couche 1
2250!----------------------------------------------------------------------------------------
2251!!! jyg le 07/02/2012
2252!!! jyg le 01/02/2017
2253       IF (iflag_split .eq. 0) THEN
2254         yt1(:) = yt(:,1)
2255         yq1(:) = yq(:,1)
2256#ifdef ISO
2257         yxt1(:,:) = yxt(:,:,1)
2258#endif
2259
2260       ELSE IF (iflag_split .ge. 1) THEN
2261#ifdef ISO
2262        call abort_physic('pbl_surface_mod 2149','isos pas encore dans iflag_split=1',1)
2263#endif
2264
2265!
2266! Cdragq computation
2267! ------------------
2268    !******************************************************************************
2269    ! Cdragq computed from cdrag
2270    ! The difference comes only from a factor (f_z0qh_oce) on z0, so that
2271    ! it can be computed inside wx_pbl0_merge
2272    ! More complicated appraches may require the propagation through
2273    ! pbl_surface of an independant cdragq variable.
2274    !******************************************************************************
2275!
2276    IF ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce) THEN
2277       ! Si on suit les formulations par exemple de Tessel, on
2278       ! a z0h=0.4*nu/u*, z0q=0.62*nu/u*, d'ou f_z0qh_oce=0.62/0.4=1.55
2279!!       ycdragq_x(1:knon)=ycdragh_x(1:knon)*                                      &
2280!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
2281!!       ycdragq_w(1:knon)=ycdragh_w(1:knon)*                                      &
2282!!            log(z1lay(1:knon)/yz0h(1:knon))/log(z1lay(1:knon)/(f_z0qh_oce*yz0h(1:knon)))
2283!
2284       DO j = 1,knon
2285         z1lay = zgeo1(j)/RG
2286         fact_cdrag = log(z1lay/yz0h(j))/log(z1lay/(f_z0qh_oce*yz0h(j)))
2287         ycdragq_x(j)=ycdragh_x(j)*fact_cdrag
2288         ycdragq_w(j)=ycdragh_w(j)*fact_cdrag
2289!!     Print *,'YYYYpbl0: fact_cdrag ', fact_cdrag
2290       ENDDO  ! j = 1,knon
2291!
2292!!  Print *,'YYYYpbl0: z1lay, yz0h, f_z0qh_oce, ycdragh_w, ycdragq_w ', &
2293!!                z1lay, yz0h(1:knon), f_z0qh_oce, ycdragh_w(1:knon), ycdragq_w(1:knon)
2294    ELSE
2295       ycdragq_x(1:knon)=ycdragh_x(1:knon)
2296       ycdragq_w(1:knon)=ycdragh_w(1:knon)
2297    ENDIF  ! ( f_z0qh_oce .ne. 1. .and. nsrf .eq.is_oce)
2298!
2299         CALL wx_pbl_prelim_0(knon, nsrf, dtime, ypplay, ypaprs, ywake_s,  &
2300                         yts, y_delta_tsurf, ygustiness, &
2301                         yt_x, yt_w, yq_x, yq_w, &
2302                         yu_x, yu_w, yv_x, yv_w, &
2303                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
2304                         ycdragm_x, ycdragm_w, &
2305                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
2306                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
2307                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2308                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2309                         Kech_h_x, Kech_h_w, Kech_h  &
2310                         )
2311         CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
2312                         BcoefQ_x, BcoefQ_w  &
2313                         )
2314         CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
2315                         ywake_s, ydTs0, ydqs0, &
2316                         yt_x, yt_w, yq_x, yq_w, &
2317                         yu_x, yu_w, yv_x, yv_w, &
2318                         ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
2319                         ycdragm_x, ycdragm_w, &
2320                         AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
2321                         AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
2322                         BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2323                         BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2324                         AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
2325                         BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
2326                         ycdragh, ycdragq, ycdragm, &
2327                         yt1, yq1, yu1, yv1 &
2328                         )
2329         IF (iflag_split .eq. 2 .AND. nsrf .ne. is_oce) THEN
2330           CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
2331                           ywake_s, ybeta, ywake_cstar, ywake_dens, &
2332                           AcoefH_x, AcoefH_w, &
2333                           BcoefH_x, BcoefH_w, &
2334                           AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
2335                           AcoefH, AcoefQ, BcoefH, BcoefQ,  &
2336                           HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2337                           phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2338                           yg_T, yg_Q, &
2339                           yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2340                           ydTs_ins, ydqs_ins &
2341                           )
2342         ELSE !
2343           AcoefH(:) = AcoefH_0(:)
2344           AcoefQ(:) = AcoefQ_0(:)
2345           BcoefH(:) = BcoefH_0(:)
2346           BcoefQ(:) = BcoefQ_0(:)
2347           yg_T(:) = 0.
2348           yg_Q(:) = 0.
2349           yGamma_dTs_phiT(:) = 0.
2350           yGamma_dQs_phiQ(:) = 0.
2351           ydTs_ins(:) = 0.
2352           ydqs_ins(:) = 0.
2353         ENDIF   ! (iflag_split .eq. 2)
2354       ENDIF  ! (iflag_split .eq.0)
2355!!!
2356       IF (prt_level >=10) THEN
2357         DO i = 1, min(1,knon)
2358           PRINT *,'pbl_surface (merge->): yt(1,:) ',yt(i,:)
2359           PRINT *,'pbl_surface (merge->): yq(1,:) ',yq(i,:)
2360           PRINT *,'pbl_surface (merge->): yu(1,:) ',yu(i,:)
2361           PRINT *,'pbl_surface (merge->): yv(1,:) ',yv(i,:)
2362           PRINT *,'pbl_surface (merge->): AcoefH(1), AcoefQ(1), AcoefU(1), AcoefV(1) ', &
2363                                           AcoefH(i), AcoefQ(i), AcoefU(i), AcoefV(i)
2364           PRINT *,'pbl_surface (merge->): BcoefH(1), BcoefQ(1), BcoefU(1), BcoefV(1) ', &
2365                                           BcoefH(i), BcoefQ(i), BcoefU(i), BcoefV(i)
2366         ENDDO
2367
2368       ENDIF
2369
2370!  Save initial value of z0h for use in evappot (z0h wiil be computed again in the surface models)
2371          yz0h_old(1:knon) = yz0h(1:knon)
2372!
2373!****************************************************************************************
2374!
2375! Calulate t2m and q2m for the case of calculation at land grid points
2376! t2m and q2m are needed as input to ORCHIDEE
2377!
2378!****************************************************************************************
2379       IF (nsrf == is_ter) THEN
2380
2381          DO i = 1, knon
2382             zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
2383                  * (ypaprs(i,1)-ypplay(i,1))
2384          ENDDO
2385
2386          ! Calculate the temperature et relative humidity at 2m and the wind at 10m
2387          IF (iflag_new_t2mq2m==1) THEN
2388           CALL stdlevvarn(klon, knon, is_ter, zxli, &
2389               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
2390               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
2391               yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
2392               yn2mout(:, nsrf, :))
2393          ELSE
2394          CALL stdlevvar(klon, knon, is_ter, zxli, &
2395               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
2396               yts, yqsurf, yz0m, yz0h, ypaprs(:,1), ypplay(:,1), &
2397               yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
2398          ENDIF
2399         
2400       ENDIF
2401
2402!****************************************************************************************
2403!
2404! 10) Switch according to current surface
2405!     It is necessary to start with the continental surfaces because the ocean
2406!     needs their run-off.
2407!
2408!****************************************************************************************
2409       SELECT CASE(nsrf)
2410     
2411       CASE(is_ter)
2412!          print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
2413          CALL surf_land(itap, dtime, date0, jour, knon, ni,&
2414               rlon, rlat, yrmu0, &
2415               debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
2416!!jyg               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2417               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
2418               AcoefH, AcoefQ, BcoefH, BcoefQ, &
2419               AcoefU, AcoefV, BcoefU, BcoefV, &
2420               ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2421               ylwdown, yq2m, yt2m, &
2422               ysnow, yqsol, yagesno, ytsoil, &
2423               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,yfluxbs,&
2424               yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
2425               y_flux_u1, y_flux_v1, &
2426               yveget,ylai,yheight   &
2427#ifdef ISO
2428         &      ,yxtrain_f, yxtsnow_f,yxt1, &
2429         &      yxtsnow,yxtsol,yxtevap,h1, &
2430         &      yrunoff_diag,yxtrunoff_diag,yRland_ice &
2431#endif               
2432         &      )
2433 
2434!FC quid qd yveget ylai yheight ne sont pas definit
2435!FC  yveget,ylai,yheight, &
2436            IF (ifl_pbltree .ge. 1) THEN
2437              CALL   freinage(knon, yu, yv, yt, &
2438!                yveget,ylai, yheight,ypaprs,ypplay,y_d_u_frein,y_d_v_frein)
2439                yveget,ylai, yheight,ypaprs,ypplay,y_treedrg, y_d_u_frein,y_d_v_frein)
2440            ENDIF
2441
2442               
2443! Special DICE MPL 05082013 puis BOMEX
2444       IF (ok_prescr_ust) THEN
2445          DO j=1,knon
2446!         ysnow(:)=0.
2447!         yqsol(:)=0.
2448!         yagesno(:)=50.
2449!         ytsoil(:,:)=300.
2450!         yz0_new(:)=0.001
2451!         yevap(:)=flat/RLVTT
2452!         yfluxlat(:)=-flat
2453!         yfluxsens(:)=-fsens
2454!         yqsurf(:)=0.
2455!         ytsurf_new(:)=tg
2456!         y_dflux_t(:)=0.
2457!         y_dflux_q(:)=0.
2458          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)
2459          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)
2460          ENDDO
2461      ENDIF
2462
2463#ifdef ISOVERIF
2464        DO j=1,knon
2465          DO ixt=1,ntraciso
2466            CALL iso_verif_noNaN(yxtevap(ixt,j), &
2467         &      'pbl_surface 1056a: apres surf_land')
2468          ENDDO
2469          DO ixt=1,niso
2470            CALL iso_verif_noNaN(yxtsol(ixt,j), &
2471         &      'pbl_surface 1056b: apres surf_land')
2472          ENDDO
2473        ENDDO
2474#endif
2475#ifdef ISOVERIF
2476!        write(*,*) 'pbl_surface_mod 1038: sortie surf_land'
2477        DO j=1,knon
2478          IF (iso_eau >= 0) THEN     
2479                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
2480     &                                  ysnow(j),'pbl_surf_mod 1043')
2481          ENDIF !if (iso_eau.gt.0) then
2482        ENDDO !DO i=1,klon
2483#endif
2484   
2485       CASE(is_lic)
2486          ! Martin
2487
2488          IF (landice_opt .LT. 2) THEN
2489             ! Land ice is treated by LMDZ and not by ORCHIDEE
2490             CALL surf_landice(itap, dtime, knon, ni, &
2491                  rlon, rlat, debut, lafin, &
2492                  yrmu0, ylwdown, yalb, zgeo1, &
2493                  ysolsw, ysollw, yts, ypplay(:,1), &
2494                  ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt1, yq1,&
2495                  AcoefH, AcoefQ, BcoefH, BcoefQ, &
2496                  AcoefU, AcoefV, BcoefU, BcoefV, &
2497                  AcoefQBS, BcoefQBS, &
2498                  ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2499                  ysnow, yqsurf, yqsol,yqbs1, yagesno, &
2500                  ytsoil, yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap,yfluxsens,yfluxlat, &
2501                  yfluxbs, ytsurf_new, y_dflux_t, y_dflux_q, &
2502                  yzmea, yzsig, ycldt, &
2503                  ysnowhgt, yqsnow, ytoice, ysissnow, &
2504                  yalb3_new, yrunoff, &
2505                  y_flux_u1, y_flux_v1 &
2506#ifdef ISO
2507                  &    ,yxtrain_f, yxtsnow_f,yxt1,yRland_ice &
2508                  &    ,yxtsnow,yxtsol,yxtevap &
2509#endif             
2510                  &    )
2511             
2512             !jyg<
2513             !!          alb3_lic(:)=0.
2514             !>jyg
2515             DO j = 1, knon
2516                i = ni(j)
2517                alb3_lic(i) = yalb3_new(j)
2518                snowhgt(i)   = ysnowhgt(j)
2519                qsnow(i)     = yqsnow(j)
2520                to_ice(i)    = ytoice(j)
2521                sissnow(i)   = ysissnow(j)
2522                runoff(i)    = yrunoff(j)
2523             ENDDO
2524             ! Martin
2525             ! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2526             IF (ok_prescr_ust) THEN
2527                DO j=1,knon
2528                   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)
2529                   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)
2530                ENDDO
2531             ENDIF
2532
2533#ifdef ISOVERIF
2534             DO j=1,knon
2535               DO ixt=1,ntraciso
2536                 CALL iso_verif_noNaN(yxtevap(ixt,j), &
2537                        &             'pbl_surface 1095a: apres surf_landice')
2538               ENDDO
2539                do ixt=1,niso
2540                   call iso_verif_noNaN(yxtsol(ixt,j), &
2541                        &      'pbl_surface 1095b: apres surf_landice')
2542                enddo
2543             enddo
2544#endif
2545#ifdef ISOVERIF
2546             !write(*,*) 'pbl_surface_mod 1060: sortie surf_landice'
2547             do j=1,knon
2548               IF (iso_eau >= 0) THEN     
2549                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
2550                        &               ysnow(j),'pbl_surf_mod 1064')
2551               ENDIF !if (iso_eau >= 0) THEN
2552             ENDDO !DO i=1,klon
2553#endif
2554           
2555          END IF
2556         
2557       CASE(is_oce)
2558           CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb_vis, &
2559               ywindsp, rmu0, yfder, yts, &
2560               itap, dtime, jour, knon, ni, &
2561!!jyg               ypplay(:,1), zgeo1/RG, ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2562               ypplay(:,1), zgeo1(1:knon)/RG, ycdragh, ycdragm, yrain_f, ysnow_f, ybs_f, yt(:,1), yq(:,1),&    ! ym missing init
2563               AcoefH, AcoefQ, BcoefH, BcoefQ, &
2564               AcoefU, AcoefV, BcoefU, BcoefV, &
2565               ypsref, yu1, yv1, ygustiness, yrugoro, pctsrf, &
2566               ysnow, yqsurf, yagesno, &
2567               yz0m, yz0h, SFRWL,yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
2568               ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
2569               y_flux_u1, y_flux_v1, ydelta_sst(:knon), ydelta_sal(:knon), &
2570               yds_ns(:knon), ydt_ns(:knon), ydter(:knon), ydser(:knon), &
2571               ydt_ds(:knon), ytkt(:knon), ytks(:knon), ytaur(:knon), ysss &
2572#ifdef ISO
2573         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
2574         &      yxtsnow,yxtevap,h1 &
2575#endif               
2576         &      )
2577      IF (prt_level >=10) THEN
2578          print *,'arg de surf_ocean: ycdragh ',ycdragh(1:knon)
2579          print *,'arg de surf_ocean: ycdragm ',ycdragm(1:knon)
2580          print *,'arg de surf_ocean: yt ', yt(1:knon,:)
2581          print *,'arg de surf_ocean: yq ', yq(1:knon,:)
2582          print *,'arg de surf_ocean: yts ', yts(1:knon)
2583          print *,'arg de surf_ocean: AcoefH ',AcoefH(1:knon)
2584          print *,'arg de surf_ocean: AcoefQ ',AcoefQ(1:knon)
2585          print *,'arg de surf_ocean: BcoefH ',BcoefH(1:knon)
2586          print *,'arg de surf_ocean: BcoefQ ',BcoefQ(1:knon)
2587          print *,'arg de surf_ocean: yevap ',yevap(1:knon)
2588          print *,'arg de surf_ocean: yfluxsens ',yfluxsens(1:knon)
2589          print *,'arg de surf_ocean: yfluxlat ',yfluxlat(1:knon)
2590          print *,'arg de surf_ocean: ytsurf_new ',ytsurf_new(1:knon)
2591       ENDIF
2592! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2593       IF (ok_prescr_ust) THEN
2594          DO j=1,knon
2595          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)
2596          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)
2597          ENDDO
2598      ENDIF
2599         
2600       CASE(is_sic)
2601          CALL surf_seaice( &
2602!albedo SB >>>
2603               rlon, rlat, ysolsw, ysollw, yalb_vis, yfder, &
2604!albedo SB <<<
2605               itap, dtime, jour, knon, ni, &
2606               lafin, &
2607!!jyg               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
2608               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt1, yq1,&
2609               AcoefH, AcoefQ, BcoefH, BcoefQ, &
2610               AcoefU, AcoefV, BcoefU, BcoefV, &
2611               ypsref, yu1, yv1, ygustiness, pctsrf, &
2612               ysnow, yqsurf, yqsol, yagesno, ytsoil, &
2613!albedo SB >>>
2614               yz0m, yz0h, SFRWL, yalb_dir_new, yalb_dif_new, yevap, yfluxsens,yfluxlat,&
2615!albedo SB <<<
2616               ytsurf_new, y_dflux_t, y_dflux_q, &
2617               y_flux_u1, y_flux_v1 &
2618#ifdef ISO
2619         &      ,yxtrain_f, yxtsnow_f,yxt1,Roce, &
2620         &      yxtsnow,yxtsol,yxtevap,Rland_ice &
2621#endif               
2622         &      )
2623         
2624! Special DICE MPL 05082013 puis BOMEX MPL 20150410
2625       IF (ok_prescr_ust) THEN
2626          DO j=1,knon
2627          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)
2628          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)
2629          ENDDO
2630       ENDIF
2631
2632#ifdef ISOVERIF
2633        DO j=1,knon
2634          DO ixt=1,ntraciso
2635            CALL iso_verif_noNaN(yxtevap(ixt,j), &
2636         &                       'pbl_surface 1165a: apres surf_seaice')
2637          ENDDO
2638          DO ixt=1,niso
2639            CALL iso_verif_noNaN(yxtsol(ixt,j), &
2640         &      'pbl_surface 1165b: apres surf_seaice')
2641          ENDDO
2642        ENDDO
2643#endif
2644#ifdef ISOVERIF
2645        !write(*,*) 'pbl_surface_mod 1077: sortie surf_seaice'
2646        DO j=1,knon
2647          IF (iso_eau >= 0) THEN     
2648                 CALL iso_verif_egalite(yxtsnow(iso_eau,j), &
2649     &                                  ysnow(j),'pbl_surf_mod 1106')
2650          ENDIF !IF (iso_eau >= 0) THEN
2651        ENDDO !DO i=1,klon
2652#endif
2653
2654       CASE DEFAULT
2655          WRITE(lunout,*) 'Surface index = ', nsrf
2656          abort_message = 'Surface index not valid'
2657          CALL abort_physic(modname,abort_message,1)
2658       END SELECT
2659
2660
2661!****************************************************************************************
2662! 11) - Calcul the increment of surface temperature
2663!
2664!****************************************************************************************
2665
2666       IF (evap0>=0.) THEN
2667          yevap(1:knon)=evap0
2668          yevap(1:knon)=RLVTT*evap0
2669       ENDIF
2670
2671       y_d_ts(1:knon)   = ytsurf_new(1:knon) - yts(1:knon)
2672 
2673!****************************************************************************************
2674!
2675! 12) "La remontee" - "The uphill"
2676!
2677!  The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
2678!  for X=H, Q, U and V, for all vertical levels.
2679!
2680!****************************************************************************************
2681!!
2682!!!
2683!!! jyg le 10/04/2013 et EV 10/2020
2684
2685        IF (ok_forc_tsurf) THEN
2686            DO j=1,knon
2687                ytsurf_new(j)=tg
2688                y_d_ts(j) = ytsurf_new(j) - yts(j)
2689            ENDDO
2690        ENDIF ! ok_forc_tsurf
2691
2692!!!
2693        IF (ok_flux_surf) THEN
2694          IF (prt_level >=10) THEN
2695           PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
2696          ENDIF
2697          y_flux_t1(:) =  fsens
2698          y_flux_q1(:) =  flat/RLVTT
2699          yfluxlat(:) =  flat
2700!
2701!!  Test sur iflag_split retire le 2/02/2018, sans vraiment comprendre la raison de ce test. (jyg)
2702!!          IF (iflag_split .eq.0) THEN
2703             DO j=1,knon
2704             Kech_h(j) = ycdragh(j) * (1.0+SQRT(yu(j,1)**2+yv(j,1)**2)) * &
2705                  ypplay(j,1)/(RD*yt(j,1))
2706             ENDDO
2707!!          ENDIF ! (iflag_split .eq.0)
2708
2709          DO j = 1, knon
2710            yt1_new=(1./RCPD)*(AcoefH(j)+BcoefH(j)*y_flux_t1(j)*dtime)
2711            ytsurf_new(j)=yt1_new-y_flux_t1(j)/(Kech_h(j)*RCPD)
2712            ! for cases forced in flux and for which forcing in Ts is needed
2713            ! to prevent the latter to reach unrealistic value (even if not used,
2714            ! Ts is calculated and hgardfou can appear during the calculation
2715            ! of surface saturation humidity for example
2716            if (ok_forc_tsurf) ytsurf_new(j)=tg
2717          ENDDO
2718
2719          DO j=1,knon
2720          y_d_ts(j) = ytsurf_new(j) - yts(j)
2721          ENDDO
2722
2723        ELSE ! (ok_flux_surf)
2724          DO j=1,knon
2725          y_flux_t1(j) =  yfluxsens(j)
2726          y_flux_q1(j) = -yevap(j)
2727#ifdef ISO
2728          y_flux_xt1(:,:) = -yxtevap(:,:)
2729#endif
2730          ENDDO
2731        ENDIF ! (ok_flux_surf)
2732
2733        ! flux of blowing snow at the first level
2734        IF (ok_bs) THEN
2735        DO j=1,knon
2736        y_flux_bs(j)=yfluxbs(j)
2737        ENDDO
2738        ENDIF
2739!
2740! ------------------------------------------------------------------------------
2741! 12a)  Splitting
2742! ------------------------------------------------------------------------------
2743
2744       IF (iflag_split .GE. 1) THEN
2745#ifdef ISO
2746        call abort_physic('pbl_surface_mod 2607','isos pas encore dans iflag_split=1',1)
2747#endif
2748!
2749!
2750         IF (nsrf .ne. is_oce) THEN
2751!
2752!         Compute potential evaporation and aridity factor  (jyg, 20200328)
2753          ybeta_prev(:) = ybeta(:)
2754             DO j = 1, knon
2755               yqa(j) = AcoefQ(j) - BcoefQ(j)*yevap(j)*dtime
2756             ENDDO
2757!
2758          CALL wx_evappot(knon, yqa, yTsurf_new, yevap_pot)
2759!
2760          ybeta(1:knon) = min(yevap(1:knon)/yevap_pot(1:knon), 1.)
2761         
2762          IF (prt_level >=10) THEN
2763           DO j=1,knon
2764            print*,'y_flux_t1,yfluxlat,wakes' &
2765 &                ,  y_flux_t1(j), yfluxlat(j), ywake_s(j)
2766            print*,'beta_prev, beta, ytsurf_new', ybeta_prev(j), ybeta(j), ytsurf_new(j)
2767            print*,'inertia,facteur,cstar', inertia, facteur,wake_cstar(j)
2768           ENDDO
2769          ENDIF  ! (prt_level >=10)
2770!
2771! Second call to wx_pbl0_merge and wx_pbl_dts_merge in order to take into account
2772! the update of the aridity coeficient beta.
2773!
2774        CALL wx_pbl_prelim_beta(knon, dtime, ywake_s, ybeta,  &
2775                        BcoefQ_x, BcoefQ_w  &
2776                        )
2777        CALL wx_pbl0_merge(knon, ypplay, ypaprs,  &
2778                          ywake_s, ydTs0, ydqs0, &
2779                          yt_x, yt_w, yq_x, yq_w, &
2780                          yu_x, yu_w, yv_x, yv_w, &
2781                          ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
2782                          ycdragm_x, ycdragm_w, &
2783                          AcoefH_x, AcoefH_w, AcoefQ_x, AcoefQ_w, &
2784                          AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
2785                          BcoefH_x, BcoefH_w, BcoefQ_x, BcoefQ_w, &
2786                          BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
2787                          AcoefH_0, AcoefQ_0, AcoefU, AcoefV, &
2788                          BcoefH_0, BcoefQ_0, BcoefU, BcoefV, &
2789                          ycdragh, ycdragq, ycdragm, &
2790                          yt1, yq1, yu1, yv1 &
2791                          )
2792          IF (iflag_split .eq. 2) THEN
2793            CALL wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
2794                            ywake_s, ybeta, ywake_cstar, ywake_dens, &
2795                            AcoefH_x, AcoefH_w, &
2796                            BcoefH_x, BcoefH_w, &
2797                            AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
2798                            AcoefH, AcoefQ, BcoefH, BcoefQ,  &
2799                            HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2800                            phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2801                            yg_T, yg_Q, &
2802                            yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2803                            ydTs_ins, ydqs_ins &
2804                            )
2805          ELSE !
2806            AcoefH(:) = AcoefH_0(:)
2807            AcoefQ(:) = AcoefQ_0(:)
2808            BcoefH(:) = BcoefH_0(:)
2809            BcoefQ(:) = BcoefQ_0(:)
2810            yg_T(:) = 0.
2811            yg_Q(:) = 0.
2812            yGamma_dTs_phiT(:) = 0.
2813            yGamma_dQs_phiQ(:) = 0.
2814            ydTs_ins(:) = 0.
2815            ydqs_ins(:) = 0.
2816          ENDIF   ! (iflag_split .eq. 2)
2817!
2818        ELSE    ! (nsrf .ne. is_oce)
2819          ybeta(1:knon) = 1.
2820          yevap_pot(1:knon) = yevap(1:knon)
2821          AcoefH(:) = AcoefH_0(:)
2822          AcoefQ(:) = AcoefQ_0(:)
2823          BcoefH(:) = BcoefH_0(:)
2824          BcoefQ(:) = BcoefQ_0(:)
2825          yg_T(:) = 0.
2826          yg_Q(:) = 0.
2827          yGamma_dTs_phiT(:) = 0.
2828          yGamma_dQs_phiQ(:) = 0.
2829          ydTs_ins(:) = 0.
2830          ydqs_ins(:) = 0.
2831        ENDIF   ! (nsrf .ne. is_oce)
2832!
2833        CALL wx_pbl_split(knon, nsrf, dtime, ywake_s, ybeta, iflag_split, &
2834                       yg_T, yg_Q, &
2835                       yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2836                       ydTs_ins, ydqs_ins, &
2837                       y_flux_t1, y_flux_q1, y_flux_u1, y_flux_v1, &
2838!!!!                       HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
2839                       phiQ0_b, phiT0_b, &
2840                       y_flux_t1_x, y_flux_t1_w, &
2841                       y_flux_q1_x, y_flux_q1_w, &
2842                       y_flux_u1_x, y_flux_u1_w, &
2843                       y_flux_v1_x, y_flux_v1_w, &
2844                       yfluxlat_x, yfluxlat_w, &
2845                       y_delta_qsats, &
2846                       y_delta_tsurf_new, y_delta_qsurf &
2847                       )
2848!
2849         CALL wx_pbl_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
2850                       yTs, y_delta_tsurf,  &
2851                       yqsurf, yTsurf_new,  &
2852                       y_delta_tsurf_new, y_delta_qsats,  &
2853                       AcoefH_x, AcoefH_w, &
2854                       BcoefH_x, BcoefH_w, &
2855                       AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
2856                       AcoefH, AcoefQ, BcoefH, BcoefQ,  &
2857                       y_flux_t1, y_flux_q1,  &
2858                       y_flux_t1_x, y_flux_t1_w, &
2859                       y_flux_q1_x, y_flux_q1_w)
2860!
2861         IF (nsrf .ne. is_oce) THEN
2862           CALL wx_pbl_dts_check(knon, dtime, ypplay, ypaprs, ywake_s, ybeta, iflag_split, &
2863                         yTs, y_delta_tsurf,  &
2864                         yqsurf, yTsurf_new,  &
2865                         y_delta_qsats, y_delta_tsurf_new, y_delta_qsurf,  &
2866                         AcoefH_x, AcoefH_w, &
2867                         BcoefH_x, BcoefH_w, &
2868                         AcoefH_0, AcoefQ_0, BcoefH_0, BcoefQ_0,  &
2869                         AcoefH, AcoefQ, BcoefH, BcoefQ,  &
2870                         HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
2871                         phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
2872                         yg_T, yg_Q, &
2873                         yGamma_dTs_phiT, yGamma_dQs_phiQ, &
2874                         ydTs_ins, ydqs_ins, &
2875                         y_flux_t1, y_flux_q1,  &
2876                         y_flux_t1_x, y_flux_t1_w, &
2877                         y_flux_q1_x, y_flux_q1_w )
2878         ENDIF   ! (nsrf .ne. is_oce)
2879!
2880       ELSE  ! (iflag_split .ge. 1)
2881         ybeta(1:knon) = 1.
2882         yevap_pot(1:knon) = yevap(1:knon)
2883       ENDIF  ! (iflag_split .ge. 1)
2884!
2885       IF (prt_level >= 10) THEN
2886         print *,'pbl_surface, ybeta , yevap, yevap_pot ', &
2887                               ybeta(1:knon) , yevap(1:knon), yevap_pot(1:knon)
2888       ENDIF  ! (prt_level >= 10)
2889!
2890!>jyg
2891!
2892 
2893!!jyg!!   A reprendre apres reflexion   ===============================================
2894!!jyg!!
2895!!jyg!!        DO j=1,knon
2896!!jyg!!!!! nrlmd le 13/06/2011
2897!!jyg!!
2898!!jyg!!!----Diffusion dans le sol dans le cas continental seulement
2899!!jyg!!       IF (nsrf.eq.is_ter) THEN
2900!!jyg!!!----Calcul du coefficient delta_coeff
2901!!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)))
2902!!jyg!!
2903!!jyg!!!          delta_coef(j)=dtime/(inertia*sqrt(tau_eq(j)))
2904!!jyg!!          delta_coef(j)=facteur*sqrt(tau_eq(j))/inertia
2905!!jyg!!!          delta_coef(j)=0.
2906!!jyg!!       ELSE
2907!!jyg!!         delta_coef(j)=0.
2908!!jyg!!       ENDIF
2909!!jyg!!
2910!!jyg!!!----Calcul de delta_tsurf
2911!!jyg!!         y_delta_tsurf(j)=delta_coef(j)*y_delta_flux_t1(j)
2912!!jyg!!
2913!!jyg!!!----Si il n'y a pas des poches...
2914!!jyg!!         IF (wake_cstar(j).le.0.01) THEN
2915!!jyg!!           y_delta_tsurf(j)=0.
2916!!jyg!!           y_delta_flux_t1(j)=0.
2917!!jyg!!         ENDIF
2918!!jyg!!
2919!!jyg!!!-----Calcul de ybeta (evap_r\'eelle/evap_potentielle)
2920!!jyg!!!!!!! jyg le 23/02/2012
2921!!jyg!!!!!!!
2922!!jyg!!!!        ybeta(j)=y_flux_q1(j)   /    &
2923!!jyg!!!! &        (Kech_h(j)*(yq(j,1)-yqsatsurf(j)))
2924!!jyg!!!!!!        ybeta(j)=-1.*yevap(j)   /    &
2925!!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)))
2926!!jyg!!!!!!! fin jyg
2927!!jyg!!!!!
2928!!jyg!!
2929!!jyg!!       ENDDO
2930!!jyg!!
2931!!jyg!!!!! fin nrlmd le 13/06/2011
2932!!jyg!!
2933       IF (iflag_split .ge. 1) THEN
2934       IF (prt_level >=10) THEN
2935        DO j = 1, knon
2936         print*,'Chx,Chw,Ch', ycdragh_x(j), ycdragh_w(j), ycdragh(j)
2937         print*,'Khx,Khw,Kh', Kech_h_x(j), Kech_h_w(j), Kech_h(j)
2938         print*,'t1x, t1w, t1, t1_ancien', &
2939 &               yt_x(j,1), yt_w(j,1),  yt(j,1), t(j,1)
2940         print*,'delta_coef,delta_flux,delta_tsurf,tau', delta_coef(j), y_delta_flux_t1(j), y_delta_tsurf(j), tau_eq(j)
2941        ENDDO
2942
2943        DO j=1,knon
2944         print*,'fluxT_x, fluxT_w, y_flux_t1, fluxQ_x, fluxQ_w, yfluxlat, wakes' &
2945 &             , 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)
2946         print*,'beta, ytsurf_new ', ybeta(j), ytsurf_new(j)
2947         print*,'inertia, facteur, cstar', inertia, facteur,wake_cstar(j)
2948        ENDDO
2949       ENDIF  ! (prt_level >=10)
2950
2951!!! jyg le 07/02/2012
2952       ENDIF  ! (iflag_split .ge.1)
2953!!!
2954
2955!!! jyg le 07/02/2012
2956       IF (iflag_split .eq.0) THEN
2957!!!
2958!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
2959        CALL climb_hq_up(knon, dtime, yt, yq, &
2960            y_flux_q1, y_flux_t1, ypaprs, ypplay, &
2961!!! jyg le 07/02/2012
2962            AcoefH, AcoefQ, BcoefH, BcoefQ, &
2963            CcoefH, CcoefQ, DcoefH, DcoefQ, &
2964            Kcoef_hq, gama_q, gama_h, &
2965!!!
2966            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:) &
2967#ifdef ISO
2968        &    ,yxt,y_flux_xt1 &
2969        &    ,AcoefXT,BcoefXT,CcoefXT,DcoefXT,gama_xt &
2970        &    ,y_flux_xt(:,:,:),y_d_xt(:,:,:) &
2971#endif
2972        &    )   
2973       ELSE  !(iflag_split .eq.0)
2974        CALL climb_hq_up(knon, dtime, yt_x, yq_x, &
2975            y_flux_q1_x, y_flux_t1_x, ypaprs, ypplay, &
2976!!! nrlmd le 02/05/2011
2977            AcoefH_x, AcoefQ_x, BcoefH_x, BcoefQ_x, &
2978            CcoefH_x, CcoefQ_x, DcoefH_x, DcoefQ_x, &
2979            Kcoef_hq_x, gama_q_x, gama_h_x, &
2980!!!
2981            y_flux_q_x(:,:), y_flux_t_x(:,:), y_d_q_x(:,:), y_d_t_x(:,:) &
2982#ifdef ISO
2983        &    ,yxt_x,y_flux_xt1_x &
2984        &    ,AcoefXT_x,BcoefXT_x,CcoefXT_x,DcoefXT_x,gama_xt_x &
2985        &    ,y_flux_xt_x(:,:,:),y_d_xt_x(:,:,:) &
2986#endif
2987        &    )   
2988!
2989       CALL climb_hq_up(knon, dtime, yt_w, yq_w, &
2990            y_flux_q1_w, y_flux_t1_w, ypaprs, ypplay, &
2991!!! nrlmd le 02/05/2011
2992            AcoefH_w, AcoefQ_w, BcoefH_w, BcoefQ_w, &
2993            CcoefH_w, CcoefQ_w, DcoefH_w, DcoefQ_w, &
2994            Kcoef_hq_w, gama_q_w, gama_h_w, &
2995!!!
2996            y_flux_q_w(:,:), y_flux_t_w(:,:), y_d_q_w(:,:), y_d_t_w(:,:) &
2997#ifdef ISO
2998        &    ,yxt_w,y_flux_xt1_w &
2999        &    ,AcoefXT_w,BcoefXT_w,CcoefXT_w,DcoefXT_w,gama_xt_w &
3000        &    ,y_flux_xt_w(:,:,:),y_d_xt_w(:,:,:) &
3001#endif
3002        &    )   
3003!!!
3004       ENDIF  ! (iflag_split .eq.0)
3005!!!
3006
3007!!! jyg le 07/02/2012
3008       IF (iflag_split .eq.0) THEN
3009!!!
3010!!! nrlmd & jyg les 02/05/2011, 13/06/2011, 05/02/2012
3011        CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
3012!!! jyg le 07/02/2012
3013            AcoefU, AcoefV, BcoefU, BcoefV, &
3014            CcoefU, CcoefV, DcoefU, DcoefV, &
3015            Kcoef_m, &
3016!!!
3017            y_flux_u, y_flux_v, y_d_u, y_d_v)
3018     y_d_t_diss(:,:)=0.
3019     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
3020        CALL yamada_c(knon,dtime,ypaprs,ypplay &
3021    &   ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
3022    &   ,iflag_pbl)
3023     ENDIF
3024!     print*,'yamada_c OK'
3025
3026       ELSE  !(iflag_split .eq.0)
3027        CALL climb_wind_up(knon, dtime, yu_x, yv_x, y_flux_u1_x, y_flux_v1_x, &
3028!!! nrlmd le 02/05/2011
3029            AcoefU_x, AcoefV_x, BcoefU_x, BcoefV_x, &
3030            CcoefU_x, CcoefV_x, DcoefU_x, DcoefV_x, &
3031            Kcoef_m_x, &
3032!!!
3033            y_flux_u_x, y_flux_v_x, y_d_u_x, y_d_v_x)
3034!
3035     y_d_t_diss_x(:,:)=0.
3036     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
3037        CALL yamada_c(knon,dtime,ypaprs,ypplay &
3038    &   ,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 &
3039        ,ycoefq_x,y_d_t_diss_x,yustar_x &
3040    &   ,iflag_pbl)
3041     ENDIF
3042!     print*,'yamada_c OK'
3043
3044        CALL climb_wind_up(knon, dtime, yu_w, yv_w, y_flux_u1_w, y_flux_v1_w, &
3045!!! nrlmd le 02/05/2011
3046            AcoefU_w, AcoefV_w, BcoefU_w, BcoefV_w, &
3047            CcoefU_w, CcoefV_w, DcoefU_w, DcoefV_w, &
3048            Kcoef_m_w, &
3049!!!
3050            y_flux_u_w, y_flux_v_w, y_d_u_w, y_d_v_w)
3051!!!
3052     y_d_t_diss_w(:,:)=0.
3053     IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
3054        CALL yamada_c(knon,dtime,ypaprs,ypplay &
3055    &   ,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 &
3056        ,ycoefq_w,y_d_t_diss_w,yustar_w &
3057    &   ,iflag_pbl)
3058     ENDIF
3059!     print*,'yamada_c OK'
3060!
3061        IF (prt_level >=10) THEN
3062         print *, 'After climbing up, lfuxlat_x, fluxlat_w ', &
3063               yfluxlat_x(1:knon), yfluxlat_w(1:knon)
3064        ENDIF
3065!
3066       ENDIF  ! (iflag_split .eq.0)
3067
3068       IF (ok_bs) THEN
3069            CALL climb_qbs_up(knon, dtime, yqbs, &
3070            y_flux_bs, ypaprs, ypplay, &
3071            AcoefQBS, BcoefQBS, &
3072            CcoefQBS, DcoefQBS, &
3073            Kcoef_qbs, gama_qbs, &
3074            y_flux_qbs(:,:), y_d_qbs(:,:))
3075       ENDIF
3076
3077!!!
3078!!
3079!!        DO j = 1, knon
3080!!          y_dflux_t(j) = y_dflux_t(j) * ypct(j)
3081!!          y_dflux_q(j) = y_dflux_q(j) * ypct(j)
3082!!        ENDDO
3083!!
3084!****************************************************************************************
3085! 13) Transform variables for output format :
3086!     - Decompress
3087!     - Multiply with pourcentage of current surface
3088!     - Cumulate in global variable
3089!
3090!****************************************************************************************
3091
3092
3093!!! jyg le 07/02/2012
3094       IF (iflag_split.EQ.0) THEN
3095!!!
3096        DO k = 1, klev
3097           DO j = 1, knon
3098             i = ni(j)
3099             y_d_t_diss(j,k)  = y_d_t_diss(j,k) * ypct(j)
3100             y_d_t(j,k)  = y_d_t(j,k) * ypct(j)
3101             y_d_q(j,k)  = y_d_q(j,k) * ypct(j)
3102             y_d_u(j,k)  = y_d_u(j,k) * ypct(j)
3103             y_d_v(j,k)  = y_d_v(j,k) * ypct(j)
3104!FC
3105             IF  (nsrf .EQ. is_ter .and. ifl_pbltree .GE. 1) THEN
3106!            if (y_d_u_frein(j,k).ne.0. ) then
3107!        print*, nsrf,'IS_TER ++', y_d_u_frein(j,k)*ypct(j),y_d_u(j,k),j,k
3108!            ENDIF
3109               y_d_u(j,k) =y_d_u(j,k) + y_d_u_frein(j,k)*ypct(j)
3110               y_d_v(j,k) =y_d_v(j,k) + y_d_v_frein(j,k)*ypct(j)
3111               treedrg(i,k,nsrf)=y_treedrg(j,k)
3112             ELSE
3113               treedrg(i,k,nsrf)=0.
3114             ENDIF
3115!FC
3116             flux_t(i,k,nsrf) = y_flux_t(j,k)
3117             flux_q(i,k,nsrf) = y_flux_q(j,k)
3118             flux_u(i,k,nsrf) = y_flux_u(j,k)
3119             flux_v(i,k,nsrf) = y_flux_v(j,k)
3120
3121#ifdef ISO
3122             DO ixt=1,ntraciso
3123                y_d_xt(ixt,j,k)  = y_d_xt(ixt,j,k) * ypct(j)
3124                flux_xt(ixt,i,k,nsrf) = y_flux_xt(ixt,j,k)
3125             ENDDO ! DO ixt=1,ntraciso
3126             h1_diag(i)=h1(j)
3127#endif
3128
3129           ENDDO
3130        ENDDO
3131
3132#ifdef ISO
3133#ifdef ISOVERIF
3134        if (iso_eau.gt.0) then
3135         call iso_verif_egalite_vect2D( &
3136                y_d_xt,y_d_q, &
3137                'pbl_surface_mod 2600',ntraciso,klon,klev)
3138        endif       
3139#endif
3140#endif
3141
3142       ELSE  !(iflag_split .eq.0)
3143
3144! Tendances hors poches
3145        DO k = 1, klev
3146          DO j = 1, knon
3147            i = ni(j)
3148            y_d_t_diss_x(j,k)  = y_d_t_diss_x(j,k) * ypct(j)
3149            y_d_t_x(j,k)  = y_d_t_x(j,k) * ypct(j)
3150            y_d_q_x(j,k)  = y_d_q_x(j,k) * ypct(j)
3151            y_d_u_x(j,k)  = y_d_u_x(j,k) * ypct(j)
3152            y_d_v_x(j,k)  = y_d_v_x(j,k) * ypct(j)
3153
3154            flux_t_x(i,k,nsrf) = y_flux_t_x(j,k)
3155            flux_q_x(i,k,nsrf) = y_flux_q_x(j,k)
3156            flux_u_x(i,k,nsrf) = y_flux_u_x(j,k)
3157            flux_v_x(i,k,nsrf) = y_flux_v_x(j,k)
3158
3159#ifdef ISO
3160            DO ixt=1,ntraciso
3161              y_d_xt_x(ixt,j,k)  = y_d_xt_x(ixt,j,k) * ypct(j)
3162              flux_xt_x(ixt,i,k,nsrf) = y_flux_xt_x(ixt,j,k)
3163            ENDDO ! DO ixt=1,ntraciso
3164#endif
3165          ENDDO
3166        ENDDO
3167
3168! Tendances dans les poches
3169        DO k = 1, klev
3170          DO j = 1, knon
3171            i = ni(j)
3172            y_d_t_diss_w(j,k)  = y_d_t_diss_w(j,k) * ypct(j)
3173            y_d_t_w(j,k)  = y_d_t_w(j,k) * ypct(j)
3174            y_d_q_w(j,k)  = y_d_q_w(j,k) * ypct(j)
3175            y_d_u_w(j,k)  = y_d_u_w(j,k) * ypct(j)
3176            y_d_v_w(j,k)  = y_d_v_w(j,k) * ypct(j)
3177
3178            flux_t_w(i,k,nsrf) = y_flux_t_w(j,k)
3179            flux_q_w(i,k,nsrf) = y_flux_q_w(j,k)
3180            flux_u_w(i,k,nsrf) = y_flux_u_w(j,k)
3181            flux_v_w(i,k,nsrf) = y_flux_v_w(j,k)
3182
3183#ifdef ISO
3184            DO ixt=1,ntraciso
3185              y_d_xt_w(ixt,j,k)  = y_d_xt_w(ixt,j,k) * ypct(j)
3186              flux_xt_w(ixt,i,k,nsrf) = y_flux_xt_w(ixt,j,k)
3187            ENDDO ! do ixt=1,ntraciso
3188#endif
3189
3190          ENDDO
3191        ENDDO
3192
3193! Flux, tendances et Tke moyenne dans la maille
3194        DO k = 1, klev
3195          DO j = 1, knon
3196            i = ni(j)
3197            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))
3198            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))
3199            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))
3200            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))
3201#ifdef ISO
3202            DO ixt=1,ntraciso
3203              flux_xt(ixt,i,k,nsrf) = flux_xt_x(ixt,i,k,nsrf)+ywake_s(j)*(flux_xt_w(ixt,i,k,nsrf)-flux_xt_x(ixt,i,k,nsrf))
3204            ENDDO ! do ixt=1,ntraciso
3205#endif
3206          ENDDO
3207        ENDDO
3208        DO j=1,knon
3209          yfluxlat(j)=yfluxlat_x(j)+ywake_s(j)*(yfluxlat_w(j)-yfluxlat_x(j))
3210        ENDDO
3211        IF (prt_level >=10) THEN
3212          print *,' nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf) ', &
3213                    nsrf, flux_t(:,1,nsrf), flux_t_x(:,1,nsrf), flux_t_w(:,1,nsrf)
3214        ENDIF
3215
3216        DO k = 1, klev
3217          DO j = 1, knon
3218            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))
3219            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))
3220            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))
3221            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))
3222            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))
3223          ENDDO
3224        ENDDO
3225
3226       ENDIF  ! (iflag_split .eq.0)
3227!!!
3228
3229       ! tendencies of blowing snow
3230       IF (ok_bs) THEN
3231           DO k = 1, klev   
3232            DO j = 1, knon
3233                i = ni(j)
3234                y_d_qbs(j,k)=y_d_qbs(j,k) * ypct(j)
3235                flux_qbs(i,k,nsrf) = y_flux_qbs(j,k)
3236            ENDDO
3237          ENDDO
3238       ENDIF
3239
3240
3241       DO j = 1, knon
3242          i = ni(j)
3243          evap(i,nsrf) = - flux_q(i,1,nsrf)                  !jyg
3244          if (ok_bs) then ; snowerosion(i,nsrf)=flux_qbs(i,1,nsrf); endif
3245          beta(i,nsrf) = ybeta(j)                             !jyg
3246          d_ts(i,nsrf) = y_d_ts(j)
3247!albedo SB >>>
3248          DO k=1,nsw
3249            alb_dir(i,k,nsrf) = yalb_dir_new(j,k)
3250            alb_dif(i,k,nsrf) = yalb_dif_new(j,k)
3251          ENDDO
3252!albedo SB <<<
3253          snow(i,nsrf) = ysnow(j) 
3254          qsurf(i,nsrf) = yqsurf(j)
3255          z0m(i,nsrf) = yz0m(j)
3256          z0h(i,nsrf) = yz0h(j)
3257          fluxlat(i,nsrf) = yfluxlat(j)
3258          agesno(i,nsrf) = yagesno(j) 
3259          cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
3260          cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
3261          dflux_t(i) = dflux_t(i) + y_dflux_t(j)*ypct(j)
3262          dflux_q(i) = dflux_q(i) + y_dflux_q(j)*ypct(j)
3263#ifdef ISO
3264        DO ixt=1,niso
3265          xtsnow(ixt,i,nsrf) = yxtsnow(ixt,j) 
3266        ENDDO
3267        DO ixt=1,ntraciso
3268          xtevap(ixt,i,nsrf) = - flux_xt(ixt,i,1,nsrf)
3269          dflux_xt(ixt,i) = dflux_xt(ixt,i) + y_dflux_xt(ixt,j)*ypct(j)
3270        ENDDO 
3271        IF (nsrf == is_lic) THEN
3272          DO ixt=1,niso
3273            Rland_ice(ixt,i) = yRland_ice(ixt,j) 
3274          ENDDO
3275        ENDIF !IF (nsrf == is_lic) THEN     
3276#ifdef ISOVERIF
3277        IF (iso_eau.gt.0) THEN 
3278          call iso_verif_egalite_choix(Rland_ice(iso_eau,i),1.0, &
3279     &         'pbl_surf_mod 1230',errmax,errmaxrel)
3280        ENDIF !if (iso_eau.gt.0) then
3281#endif       
3282#endif
3283       ENDDO
3284
3285!      print*,'Dans pbl OK2'
3286
3287!!! jyg le 07/02/2012
3288       IF (iflag_split .ge.1) THEN
3289!!!
3290!!! nrlmd le 02/05/2011
3291        DO j = 1, knon
3292          i = ni(j)
3293          fluxlat_x(i,nsrf) = yfluxlat_x(j)
3294          fluxlat_w(i,nsrf) = yfluxlat_w(j)
3295!!!
3296!!! nrlmd le 13/06/2011
3297!!jyg20170131          delta_tsurf(i,nsrf)=y_delta_tsurf(j)*ypct(j)
3298!!jyg20210118          delta_tsurf(i,nsrf)=y_delta_tsurf(j)
3299          delta_tsurf(i,nsrf)=y_delta_tsurf_new(j)
3300!
3301          delta_qsurf(i) = delta_qsurf(i) + y_delta_qsurf(j)*ypct(j)
3302!
3303          cdragh_x(i) = cdragh_x(i) + ycdragh_x(j)*ypct(j)
3304          cdragh_w(i) = cdragh_w(i) + ycdragh_w(j)*ypct(j)
3305          cdragm_x(i) = cdragm_x(i) + ycdragm_x(j)*ypct(j)
3306          cdragm_w(i) = cdragm_w(i) + ycdragm_w(j)*ypct(j)
3307          kh(i) = kh(i) + Kech_h(j)*ypct(j)
3308          kh_x(i) = kh_x(i) + Kech_h_x(j)*ypct(j)
3309          kh_w(i) = kh_w(i) + Kech_h_w(j)*ypct(j)
3310!!!
3311        ENDDO
3312!!!     
3313       ENDIF  ! (iflag_split .ge.1)
3314!!!
3315!!! nrlmd le 02/05/2011
3316!!jyg le 20/02/2011
3317!!        tke_x(:,:,nsrf)=0.
3318!!        tke_w(:,:,nsrf)=0.
3319!!jyg le 20/02/2011
3320!!        DO k = 1, klev+1
3321!!          DO j = 1, knon
3322!!            i = ni(j)
3323!!            wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
3324!!            tke(i,k,nsrf)   = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
3325!!          ENDDO
3326!!        ENDDO
3327!!jyg le 20/02/2011
3328!!        DO k = 1, klev+1
3329!!          DO j = 1, knon
3330!!            i = ni(j)
3331!!            tke(i,k,nsrf)=(1.-ywake_s(j))*tke_x(i,k,nsrf)+ywake_s(j)*tke_w(i,k,nsrf)
3332!!          ENDDO
3333!!        ENDDO
3334!!!
3335       IF (iflag_split .eq.0) THEN
3336        wake_dltke(:,:,nsrf) = 0.
3337        DO k = 1, klev+1
3338           DO j = 1, knon
3339              i = ni(j)
3340!jyg<
3341!!              tke(i,k,nsrf)    = ytke(j,k)
3342!!              tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
3343              tke_x(i,k,nsrf)    = ytke(j,k)
3344              tke_x(i,k,is_ave)  = tke_x(i,k,is_ave) + ytke(j,k)*ypct(j)
3345              eps_x(i,k,nsrf)    = yeps(j,k)
3346              eps_x(i,k,is_ave)  = eps_x(i,k,is_ave) + yeps(j,k)*ypct(j)
3347!>jyg
3348           ENDDO
3349        ENDDO
3350
3351       ELSE  ! (iflag_split .eq.0)
3352        DO k = 1, klev+1
3353          DO j = 1, knon
3354            i = ni(j)
3355            wake_dltke(i,k,nsrf) = ytke_w(j,k) - ytke_x(j,k)
3356!jyg<
3357!!            tke(i,k,nsrf)   = ytke_x(j,k) + ywake_s(j)*wake_dltke(i,k,nsrf)
3358!!            tke(i,k,is_ave) = tke(i,k,is_ave) + tke(i,k,nsrf)*ypct(j)
3359            tke_x(i,k,nsrf)   = ytke_x(j,k)
3360            tke_x(i,k,is_ave)   = tke_x(i,k,is_ave) + tke_x(i,k,nsrf)*ypct(j)       
3361            eps_x(i,k,nsrf)   = yeps_x(j,k)
3362            eps_x(i,k,is_ave)   = eps_x(i,k,is_ave) + eps_x(i,k,nsrf)*ypct(j)
3363            wake_dltke(i,k,is_ave)   = wake_dltke(i,k,is_ave) + wake_dltke(i,k,nsrf)*ypct(j)
3364           
3365
3366!>jyg
3367          ENDDO
3368        ENDDO
3369       ENDIF  ! (iflag_split .eq.0)
3370!!!
3371       DO k = 2, klev
3372          DO j = 1, knon
3373             i = ni(j)
3374             zcoefh(i,k,nsrf) = ycoefh(j,k)
3375             zcoefm(i,k,nsrf) = ycoefm(j,k)
3376             zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
3377             zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
3378          ENDDO
3379       ENDDO
3380
3381!      print*,'Dans pbl OK3'
3382
3383       IF ( nsrf .EQ. is_ter ) THEN
3384          DO j = 1, knon
3385             i = ni(j)
3386             qsol(i) = yqsol(j)
3387#ifdef ISO
3388             runoff_diag(i)=yrunoff_diag(j)   
3389             DO ixt=1,niso
3390               xtsol(ixt,i) = yxtsol(ixt,j)
3391               xtrunoff_diag(ixt,i)=yxtrunoff_diag(ixt,j)
3392             ENDDO
3393#endif
3394          ENDDO
3395       ENDIF
3396       
3397!jyg<
3398!!       ftsoil(:,:,nsrf) = 0.
3399!>jyg
3400       DO k = 1, nsoilmx
3401          DO j = 1, knon
3402             i = ni(j)
3403             ftsoil(i, k, nsrf) = ytsoil(j,k)
3404          ENDDO
3405       ENDDO
3406
3407#ifdef ISO
3408#ifdef ISOVERIF
3409       !write(*,*) 'pbl_surface 2858'
3410       DO i = 1, klon
3411         DO ixt=1,niso
3412           call iso_verif_noNaN(xtsol(ixt,i),'pbl_surface 1405')
3413         ENDDO
3414       ENDDO
3415#endif
3416#ifdef ISOVERIF
3417     IF (iso_eau.gt.0) THEN
3418        call iso_verif_egalite_vect2D( &
3419                y_d_xt,y_d_q, &
3420                'pbl_surface_mod 1261',ntraciso,klon,klev)
3421     ENDIF !if (iso_eau.gt.0) then
3422#endif
3423#endif
3424!!! jyg le 07/02/2012
3425       IF (iflag_split .ge.1) THEN
3426!!!
3427!!! nrlmd+jyg le 02/05/2011 et le 20/02/2012
3428        DO k = 1, klev
3429          DO j = 1, knon
3430           i = ni(j)
3431           d_t_diss_x(i,k) = d_t_diss_x(i,k) + y_d_t_diss_x(j,k)
3432           d_t_x(i,k) = d_t_x(i,k) + y_d_t_x(j,k)
3433           d_q_x(i,k) = d_q_x(i,k) + y_d_q_x(j,k)
3434           d_u_x(i,k) = d_u_x(i,k) + y_d_u_x(j,k)
3435           d_v_x(i,k) = d_v_x(i,k) + y_d_v_x(j,k)
3436!
3437           d_t_diss_w(i,k) = d_t_diss_w(i,k) + y_d_t_diss_w(j,k)
3438           d_t_w(i,k) = d_t_w(i,k) + y_d_t_w(j,k)
3439           d_q_w(i,k) = d_q_w(i,k) + y_d_q_w(j,k)
3440           d_u_w(i,k) = d_u_w(i,k) + y_d_u_w(j,k)
3441           d_v_w(i,k) = d_v_w(i,k) + y_d_v_w(j,k)
3442#ifdef ISO
3443           DO ixt=1,ntraciso
3444             d_xt_x(ixt,i,k) = d_xt_x(ixt,i,k) + y_d_xt_x(ixt,j,k)
3445             d_xt_w(ixt,i,k) = d_xt_w(ixt,i,k) + y_d_xt_w(ixt,j,k)
3446           ENDDO ! DO ixt=1,ntraciso
3447#endif
3448
3449!
3450!!           d_wake_dlt(i,k) = d_wake_dlt(i,k) + y_d_t_w(i,k)-y_d_t_x(i,k)
3451!!           d_wake_dlq(i,k) = d_wake_dlq(i,k) + y_d_q_w(i,k)-y_d_q_x(i,k)
3452          ENDDO
3453        ENDDO
3454!!!
3455       ENDIF  ! (iflag_split .ge.1)
3456!!!
3457       
3458       DO k = 1, klev
3459          DO j = 1, knon
3460             i = ni(j)
3461             d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
3462             d_t(i,k) = d_t(i,k) + y_d_t(j,k)
3463             d_q(i,k) = d_q(i,k) + y_d_q(j,k)
3464#ifdef ISO
3465             DO ixt=1,ntraciso
3466               d_xt(ixt,i,k) = d_xt(ixt,i,k) + y_d_xt(ixt,j,k)
3467             ENDDO !DO ixt=1,ntraciso
3468#endif
3469             d_u(i,k) = d_u(i,k) + y_d_u(j,k)
3470             d_v(i,k) = d_v(i,k) + y_d_v(j,k)
3471          ENDDO
3472       ENDDO
3473
3474
3475       IF (ok_bs) THEN
3476         DO k = 1, klev
3477         DO j = 1, knon
3478         i = ni(j)
3479         d_qbs(i,k) = d_qbs(i,k) + y_d_qbs(j,k)
3480         ENDDO
3481         ENDDO
3482        ENDIF
3483
3484#ifdef ISO
3485#ifdef ISOVERIF
3486!        write(*,*) 'd_q,d_xt(iso_eau,554,19)=',d_q(554,19),d_xt(iso_eau,554,19)
3487!        write(*,*) 'pbl_surface 2929: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)
3488!        write(*,*) 'y_d_q,y_d_xt(iso_eau,2,1)=',y_d_q(2,1),y_d_xt(iso_eau,2,1)
3489!        write(*,*) 'iso_eau.gt.0=',iso_eau.gt.0
3490        call iso_verif_noNaN_vect2D( &
3491     &           d_xt, &
3492     &           'pbl_surface 1385',ntraciso,klon,klev) 
3493     IF (iso_eau >= 0) THEN
3494        call iso_verif_egalite_vect2D( &
3495                y_d_xt,y_d_q, &
3496                'pbl_surface_mod 2945',ntraciso,klon,klev)
3497        call iso_verif_egalite_vect2D( &
3498                d_xt,d_q, &
3499                'pbl_surface_mod 1276',ntraciso,klon,klev)
3500     ENDIF !IF (iso_eau >= 0) THEN
3501#endif
3502#endif
3503
3504!      print*,'Dans pbl OK4'
3505
3506       IF (prt_level >=10) THEN
3507         print *, 'pbl_surface tendencies for w: d_t_w, d_t_x, d_t ', &
3508          d_t_w(1:knon,1), d_t_x(1:knon,1), d_t(1:knon,1)
3509       ENDIF
3510
3511       if (nsrf == is_oce .and. activate_ocean_skin >= 1) then
3512          delta_sal = missing_val
3513          ds_ns = missing_val
3514          dt_ns = missing_val
3515          delta_sst = missing_val
3516          dter = missing_val
3517          dser = missing_val
3518          tkt = missing_val
3519          tks = missing_val
3520          taur = missing_val
3521          sss = missing_val
3522         
3523          delta_sal(ni(:knon)) = ydelta_sal(:knon)
3524          ds_ns(ni(:knon)) = yds_ns(:knon)
3525          dt_ns(ni(:knon)) = ydt_ns(:knon)
3526          delta_sst(ni(:knon)) = ydelta_sst(:knon)
3527          dter(ni(:knon)) = ydter(:knon)
3528          dser(ni(:knon)) = ydser(:knon)
3529          tkt(ni(:knon)) = ytkt(:knon)
3530          tks(ni(:knon)) = ytks(:knon)
3531          taur(ni(:knon)) = ytaur(:knon)
3532          sss(ni(:knon)) = ysss(:knon)
3533
3534          if (activate_ocean_skin == 2 .and. type_ocean == "couple") then
3535             dt_ds = missing_val
3536             dt_ds(ni(:knon)) = ydt_ds(:knon)
3537          end if
3538       end if
3539
3540
3541!****************************************************************************************
3542! 14) Calculate the temperature and relative humidity at 2m and the wind at 10m
3543!     Call HBTM
3544!
3545!****************************************************************************************
3546!!!
3547!
3548#undef T2m     
3549#define T2m     
3550#ifdef T2m
3551! Calculations of diagnostic t,q at 2m and u, v at 10m
3552
3553!      print*,'Dans pbl OK41'
3554!      print*,'tair1,yt(:,1),y_d_t(:,1)'
3555!      print*, tair1,yt(:,1),y_d_t(:,1)
3556!!! jyg le 07/02/2012
3557       IF (iflag_split .eq.0) THEN
3558        DO j=1, knon
3559          uzon(j) = yu(j,1) + y_d_u(j,1)
3560          vmer(j) = yv(j,1) + y_d_v(j,1)
3561          tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
3562          qair1(j) = yq(j,1) + y_d_q(j,1)
3563          zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3564               * (ypaprs(j,1)-ypplay(j,1))
3565          tairsol(j) = yts(j) + y_d_ts(j)
3566          qairsol(j) = yqsurf(j)
3567        ENDDO
3568       ELSE  ! (iflag_split .eq.0)
3569        DO j=1, knon
3570          uzon_x(j) = yu_x(j,1) + y_d_u_x(j,1)
3571          vmer_x(j) = yv_x(j,1) + y_d_v_x(j,1)
3572          tair1_x(j) = yt_x(j,1) + y_d_t_x(j,1) + y_d_t_diss_x(j,1)
3573          qair1_x(j) = yq_x(j,1) + y_d_q_x(j,1)
3574          zgeo1_x(j) = RD * tair1_x(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3575               * (ypaprs(j,1)-ypplay(j,1))
3576          tairsol(j) = yts(j) + y_d_ts(j)
3577!!          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf(j)
3578          tairsol_x(j) = tairsol(j) - ywake_s(j)*y_delta_tsurf_new(j)
3579          qairsol(j) = yqsurf(j)
3580        ENDDO
3581        DO j=1, knon
3582          uzon_w(j) = yu_w(j,1) + y_d_u_w(j,1)
3583          vmer_w(j) = yv_w(j,1) + y_d_v_w(j,1)
3584          tair1_w(j) = yt_w(j,1) + y_d_t_w(j,1) + y_d_t_diss_w(j,1)
3585          qair1_w(j) = yq_w(j,1) + y_d_q_w(j,1)
3586          zgeo1_w(j) = RD * tair1_w(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
3587               * (ypaprs(j,1)-ypplay(j,1))
3588          tairsol_w(j) = tairsol(j) + (1.- ywake_s(j))*y_delta_tsurf(j)
3589          qairsol(j) = yqsurf(j)
3590        ENDDO
3591!!!     
3592       ENDIF  ! (iflag_split .eq.0)
3593!!!
3594       DO j=1, knon
3595!         i = ni(j)
3596!         yz0h_oupas(j) = yz0m(j)
3597!         IF(nsrf.EQ.is_oce) THEN
3598!            yz0h_oupas(j) = z0m(i,nsrf)
3599!         ENDIF
3600          psfce(j)=ypaprs(j,1)
3601          patm(j)=ypplay(j,1)
3602       ENDDO
3603
3604       IF (iflag_pbl_surface_t2m_bug==1) THEN
3605          yz0h_oupas(1:knon)=yz0m(1:knon)
3606       ELSE
3607          yz0h_oupas(1:knon)=yz0h(1:knon)
3608       ENDIF
3609       
3610!      print*,'Dans pbl OK42A'
3611!      print*,'tair1,yt(:,1),y_d_t(:,1)'
3612!      print*, tair1,yt(:,1),y_d_t(:,1)
3613
3614! Calculate the temperature and relative humidity at 2m and the wind at 10m
3615!!! jyg le 07/02/2012
3616       IF (iflag_split .eq.0) THEN
3617        IF (iflag_new_t2mq2m==1) THEN
3618           CALL stdlevvarn(klon, knon, nsrf, zxli, &
3619            uzon, vmer, tair1, qair1, zgeo1, &
3620            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3621            yt2m, yq2m, yt10m, yq10m, yu10m, yustar, &
3622            yn2mout(:, nsrf, :))
3623        ELSE
3624        CALL stdlevvar(klon, knon, nsrf, zxli, &
3625            uzon, vmer, tair1, qair1, zgeo1, &
3626            tairsol, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3627            yt2m, yq2m, yt10m, yq10m, yu10m, yustar, ypblh, rain_f, zxtsol)
3628        ENDIF
3629       ELSE  !(iflag_split .eq.0)
3630        IF (iflag_new_t2mq2m==1) THEN
3631         CALL stdlevvarn(klon, knon, nsrf, zxli, &
3632            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
3633            tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3634            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, &
3635            yn2mout_x(:, nsrf, :))
3636         CALL stdlevvarn(klon, knon, nsrf, zxli, &
3637            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
3638            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3639            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, &
3640            yn2mout_w(:, nsrf, :))
3641        ELSE
3642        CALL stdlevvar(klon, knon, nsrf, zxli, &
3643            uzon_x, vmer_x, tair1_x, qair1_x, zgeo1_x, &
3644            tairsol_x, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3645            yt2m_x, yq2m_x, yt10m_x, yq10m_x, yu10m_x, yustar_x, ypblh_x, rain_f, zxtsol)
3646        CALL stdlevvar(klon, knon, nsrf, zxli, &
3647            uzon_w, vmer_w, tair1_w, qair1_w, zgeo1_w, &
3648            tairsol_w, qairsol, yz0m, yz0h_oupas, psfce, patm, &
3649            yt2m_w, yq2m_w, yt10m_w, yq10m_w, yu10m_w, yustar_w, ypblh_w, rain_f, zxtsol)
3650        ENDIF
3651!!!
3652       ENDIF  ! (iflag_split .eq.0)
3653!!!
3654!!! jyg le 07/02/2012
3655       IF (iflag_split .eq.0) THEN
3656        DO j=1, knon
3657          i = ni(j)
3658          t2m(i,nsrf)=yt2m(j)
3659          q2m(i,nsrf)=yq2m(j)
3660     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3661          ustar(i,nsrf)=yustar(j)
3662          u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
3663          v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
3664!
3665          DO k = 1, 6
3666           n2mout(i,nsrf,k) = yn2mout(j,nsrf,k)
3667          END DO 
3668!
3669        ENDDO
3670       ELSE  !(iflag_split .eq.0)
3671        DO j=1, knon
3672          i = ni(j)
3673          t2m_x(i,nsrf)=yt2m_x(j)
3674          q2m_x(i,nsrf)=yq2m_x(j)
3675     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3676          ustar_x(i,nsrf)=yustar_x(j)
3677          u10m_x(i,nsrf)=(yu10m_x(j) * uzon_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
3678          v10m_x(i,nsrf)=(yu10m_x(j) * vmer_x(j))/SQRT(uzon_x(j)**2+vmer_x(j)**2)
3679!
3680          DO k = 1, 6
3681           n2mout_x(i,nsrf,k) = yn2mout_x(j,nsrf,k)
3682          END DO 
3683!
3684        ENDDO
3685        DO j=1, knon
3686          i = ni(j)
3687          t2m_w(i,nsrf)=yt2m_w(j)
3688          q2m_w(i,nsrf)=yq2m_w(j)
3689     ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
3690          ustar_w(i,nsrf)=yustar_w(j)
3691          u10m_w(i,nsrf)=(yu10m_w(j) * uzon_w(j))/SQRT(uzon_w(j)**2+vmer_w(j)**2)
3692          v10m_w(i,nsrf)=(yu10m_w(j) * vmer_w(j))/SQRT(uzon_w(j)**2+vmer_w(j)**2)
3693!
3694          ustar(i,nsrf) = ustar_x(i,nsrf) + wake_s(i)*(ustar_w(i,nsrf)-ustar_x(i,nsrf))
3695          u10m(i,nsrf) = u10m_x(i,nsrf) + wake_s(i)*(u10m_w(i,nsrf)-u10m_x(i,nsrf))
3696          v10m(i,nsrf) = v10m_x(i,nsrf) + wake_s(i)*(v10m_w(i,nsrf)-v10m_x(i,nsrf))
3697!
3698          DO k = 1, 6
3699           n2mout_w(i,nsrf,k) = yn2mout_w(j,nsrf,k)
3700          END DO 
3701!
3702        ENDDO
3703!!!
3704       ENDIF  ! (iflag_split .eq.0)
3705!!!
3706
3707!      print*,'Dans pbl OK43'
3708!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
3709!IM Ajoute dependance type surface
3710       IF (thermcep) THEN
3711!!! jyg le 07/02/2012
3712       IF (iflag_split .eq.0) THEN
3713          DO j = 1, knon
3714             i=ni(j)
3715             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
3716             zx_qs1  = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
3717             zx_qs1  = MIN(0.5,zx_qs1)
3718             zcor1   = 1./(1.-RETV*zx_qs1)
3719             zx_qs1  = zx_qs1*zcor1
3720             
3721             rh2m(i)   = rh2m(i)   + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
3722             qsat2m(i) = qsat2m(i) + zx_qs1  * pctsrf(i,nsrf)
3723          ENDDO
3724       ELSE  ! (iflag_split .eq.0)
3725          DO j = 1, knon
3726             i=ni(j)
3727             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_x(j) ))
3728             zx_qs1  = r2es * FOEEW(yt2m_x(j),zdelta1)/paprs(i,1)
3729             zx_qs1  = MIN(0.5,zx_qs1)
3730             zcor1   = 1./(1.-RETV*zx_qs1)
3731             zx_qs1  = zx_qs1*zcor1
3732             
3733             rh2m_x(i)   = rh2m_x(i)   + yq2m_x(j)/zx_qs1 * pctsrf(i,nsrf)
3734             qsat2m_x(i) = qsat2m_x(i) + zx_qs1  * pctsrf(i,nsrf)
3735          ENDDO
3736          DO j = 1, knon
3737             i=ni(j)
3738             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m_w(j) ))
3739             zx_qs1  = r2es * FOEEW(yt2m_w(j),zdelta1)/paprs(i,1)
3740             zx_qs1  = MIN(0.5,zx_qs1)
3741             zcor1   = 1./(1.-RETV*zx_qs1)
3742             zx_qs1  = zx_qs1*zcor1
3743             
3744             rh2m_w(i)   = rh2m_w(i)   + yq2m_w(j)/zx_qs1 * pctsrf(i,nsrf)
3745             qsat2m_w(i) = qsat2m_w(i) + zx_qs1  * pctsrf(i,nsrf)
3746          ENDDO
3747!!!     
3748       ENDIF  ! (iflag_split .eq.0)
3749!!!
3750       ENDIF
3751!
3752       IF (prt_level >=10) THEN
3753         print *, 'T2m, q2m, RH2m ', &
3754          t2m(1:knon,:), q2m(1:knon,:), rh2m(1:knon)
3755       ENDIF
3756
3757!   print*,'OK pbl 5'
3758!
3759!!! jyg le 07/02/2012
3760       IF (iflag_split .eq.0) THEN
3761        CALL hbtm(knon, ypaprs, ypplay, &
3762            yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
3763            y_flux_t,y_flux_q,yu,yv,yt,yq, &
3764            ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
3765            ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
3766          IF (prt_level >=10) THEN
3767       print *,' Arg. de HBTM: yt2m ',yt2m(1:knon)
3768       print *,' Arg. de HBTM: yt10m ',yt10m(1:knon)
3769       print *,' Arg. de HBTM: yq2m ',yq2m(1:knon)
3770       print *,' Arg. de HBTM: yq10m ',yq10m(1:knon)
3771       print *,' Arg. de HBTM: yustar ',yustar(1:knon)
3772       print *,' Arg. de HBTM: y_flux_t ',y_flux_t(1:knon,:)
3773       print *,' Arg. de HBTM: y_flux_q ',y_flux_q(1:knon,:)
3774       print *,' Arg. de HBTM: yu ',yu(1:knon,:)
3775       print *,' Arg. de HBTM: yv ',yv(1:knon,:)
3776       print *,' Arg. de HBTM: yt ',yt(1:knon,:)
3777       print *,' Arg. de HBTM: yq ',yq(1:knon,:)
3778          ENDIF
3779       ELSE  ! (iflag_split .eq.0)
3780        CALL HBTM(knon, ypaprs, ypplay, &
3781            yt2m_x,yt10m_x,yq2m_x,yq10m_x,yustar_x,ywstar_x, &
3782            y_flux_t_x,y_flux_q_x,yu_x,yv_x,yt_x,yq_x, &
3783            ypblh_x,ycapCL_x,yoliqCL_x,ycteiCL_x,ypblT_x, &
3784            ytherm_x,ytrmb1_x,ytrmb2_x,ytrmb3_x,ylcl_x)
3785          IF (prt_level >=10) THEN
3786       print *,' Arg. de HBTM: yt2m_x ',yt2m_x(1:knon)
3787       print *,' Arg. de HBTM: yt10m_x ',yt10m_x(1:knon)
3788       print *,' Arg. de HBTM: yq2m_x ',yq2m_x(1:knon)
3789       print *,' Arg. de HBTM: yq10m_x ',yq10m_x(1:knon)
3790       print *,' Arg. de HBTM: yustar_x ',yustar_x(1:knon)
3791       print *,' Arg. de HBTM: y_flux_t_x ',y_flux_t_x(1:knon,:)
3792       print *,' Arg. de HBTM: y_flux_q_x ',y_flux_q_x(1:knon,:)
3793       print *,' Arg. de HBTM: yu_x ',yu_x(1:knon,:)
3794       print *,' Arg. de HBTM: yv_x ',yv_x(1:knon,:)
3795       print *,' Arg. de HBTM: yt_x ',yt_x(1:knon,:)
3796       print *,' Arg. de HBTM: yq_x ',yq_x(1:knon,:)
3797          ENDIF
3798        CALL HBTM(knon, ypaprs, ypplay, &
3799            yt2m_w,yt10m_w,yq2m_w,yq10m_w,yustar_w,ywstar_w, &
3800            y_flux_t_w,y_flux_q_w,yu_w,yv_w,yt_w,yq_w, &
3801            ypblh_w,ycapCL_w,yoliqCL_w,ycteiCL_w,ypblT_w, &
3802            ytherm_w,ytrmb1_w,ytrmb2_w,ytrmb3_w,ylcl_w)
3803!!!     
3804       ENDIF  ! (iflag_split .eq.0)
3805!!!
3806       
3807!!! jyg le 07/02/2012
3808       IF (iflag_split .eq.0) THEN
3809!!!
3810        DO j=1, knon
3811          i = ni(j)
3812          pblh(i,nsrf)   = ypblh(j)
3813          wstar(i,nsrf)  = ywstar(j)
3814          plcl(i,nsrf)   = ylcl(j)
3815          capCL(i,nsrf)  = ycapCL(j)
3816          oliqCL(i,nsrf) = yoliqCL(j)
3817          cteiCL(i,nsrf) = ycteiCL(j)
3818          pblT(i,nsrf)   = ypblT(j)
3819          therm(i,nsrf)  = ytherm(j)
3820          trmb1(i,nsrf)  = ytrmb1(j)
3821          trmb2(i,nsrf)  = ytrmb2(j)
3822          trmb3(i,nsrf)  = ytrmb3(j)
3823        ENDDO
3824        IF (prt_level >=10) THEN
3825          print *, 'After HBTM: pblh ', pblh(1:knon,:)
3826          print *, 'After HBTM: plcl ', plcl(1:knon,:)
3827          print *, 'After HBTM: cteiCL ', cteiCL(1:knon,:)
3828        ENDIF
3829       ELSE  !(iflag_split .eq.0)
3830        DO j=1, knon
3831          i = ni(j)
3832          pblh_x(i,nsrf)   = ypblh_x(j)
3833          wstar_x(i,nsrf)  = ywstar_x(j)
3834          plcl_x(i,nsrf)   = ylcl_x(j)
3835          capCL_x(i,nsrf)  = ycapCL_x(j)
3836          oliqCL_x(i,nsrf) = yoliqCL_x(j)
3837          cteiCL_x(i,nsrf) = ycteiCL_x(j)
3838          pblT_x(i,nsrf)   = ypblT_x(j)
3839          therm_x(i,nsrf)  = ytherm_x(j)
3840          trmb1_x(i,nsrf)  = ytrmb1_x(j)
3841          trmb2_x(i,nsrf)  = ytrmb2_x(j)
3842          trmb3_x(i,nsrf)  = ytrmb3_x(j)
3843        ENDDO
3844        IF (prt_level >=10) THEN
3845          print *, 'After HBTM: pblh_x ', pblh_x(1:knon,:)
3846          print *, 'After HBTM: plcl_x ', plcl_x(1:knon,:)
3847          print *, 'After HBTM: cteiCL_x ', cteiCL_x(1:knon,:)
3848        ENDIF
3849        DO j=1, knon
3850          i = ni(j)
3851          pblh_w(i,nsrf)   = ypblh_w(j)
3852          wstar_w(i,nsrf)  = ywstar_w(j)
3853          plcl_w(i,nsrf)   = ylcl_w(j)
3854          capCL_w(i,nsrf)  = ycapCL_w(j)
3855          oliqCL_w(i,nsrf) = yoliqCL_w(j)
3856          cteiCL_w(i,nsrf) = ycteiCL_w(j)
3857          pblT_w(i,nsrf)   = ypblT_w(j)
3858          therm_w(i,nsrf)  = ytherm_w(j)
3859          trmb1_w(i,nsrf)  = ytrmb1_w(j)
3860          trmb2_w(i,nsrf)  = ytrmb2_w(j)
3861          trmb3_w(i,nsrf)  = ytrmb3_w(j)
3862        ENDDO
3863        IF (prt_level >=10) THEN
3864          print *, 'After HBTM: pblh_w ', pblh_w(1:knon,:)
3865          print *, 'After HBTM: plcl_w ', plcl_w(1:knon,:)
3866          print *, 'After HBTM: cteiCL_w ', cteiCL_w(1:knon,:)
3867        ENDIF
3868!!!
3869       ENDIF  ! (iflag_split .eq.0)
3870!!!
3871
3872!   print*,'OK pbl 6'
3873#else
3874! T2m not defined
3875! No calculation
3876       PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
3877#endif
3878
3879!****************************************************************************************
3880! 15) End of loop over different surfaces
3881!
3882!****************************************************************************************
3883    ENDDO loop_nbsrf
3884!
3885!----------------------------------------------------------------------------------------
3886!   Reset iflag_split
3887!
3888   iflag_split=iflag_split_ref
3889
3890#ifdef ISO
3891#ifdef ISOVERIF
3892!        write(*,*) 'pbl_surface tmp 3249: d_q,d_xt(iso_eau,2,1)=',d_q(2,1),d_xt(iso_eau,2,1)
3893    IF (iso_eau >= 0) THEN
3894        call iso_verif_egalite_vect2D( &
3895                d_xt,d_q, &
3896                'pbl_surface_mod 1276',ntraciso,klon,klev)
3897    ENDIF !IF (iso_eau >= 0) THEN
3898#endif
3899#endif
3900
3901!****************************************************************************************
3902! 16) Calculate the mean value over all sub-surfaces for some variables
3903!
3904!****************************************************************************************
3905   
3906    z0m(:,nbsrf+1) = 0.0
3907    z0h(:,nbsrf+1) = 0.0
3908    DO nsrf = 1, nbsrf
3909       DO i = 1, klon
3910          z0m(i,nbsrf+1) = z0m(i,nbsrf+1) + z0m(i,nsrf)*pctsrf(i,nsrf)
3911          z0h(i,nbsrf+1) = z0h(i,nbsrf+1) + z0h(i,nsrf)*pctsrf(i,nsrf)
3912       ENDDO
3913    ENDDO
3914
3915!   print*,'OK pbl 7'
3916    zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
3917    zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
3918    zxfluxt_x(:,:) = 0.0 ; zxfluxq_x(:,:) = 0.0
3919    zxfluxu_x(:,:) = 0.0 ; zxfluxv_x(:,:) = 0.0
3920    zxfluxt_w(:,:) = 0.0 ; zxfluxq_w(:,:) = 0.0
3921    zxfluxu_w(:,:) = 0.0 ; zxfluxv_w(:,:) = 0.0
3922#ifdef ISO
3923      zxfluxxt(:,:,:) = 0.0
3924      zxfluxxt_x(:,:,:) = 0.0
3925      zxfluxxt_w(:,:,:) = 0.0
3926#endif
3927
3928
3929!!! jyg le 07/02/2012
3930       IF (iflag_split .ge.1) THEN
3931!!!
3932!!! nrlmd & jyg les 02/05/2011, 05/02/2012
3933
3934        DO nsrf = 1, nbsrf
3935          DO k = 1, klev
3936            DO i = 1, klon
3937              zxfluxt_x(i,k) = zxfluxt_x(i,k) + flux_t_x(i,k,nsrf) * pctsrf(i,nsrf)
3938              zxfluxq_x(i,k) = zxfluxq_x(i,k) + flux_q_x(i,k,nsrf) * pctsrf(i,nsrf)
3939              zxfluxu_x(i,k) = zxfluxu_x(i,k) + flux_u_x(i,k,nsrf) * pctsrf(i,nsrf)
3940              zxfluxv_x(i,k) = zxfluxv_x(i,k) + flux_v_x(i,k,nsrf) * pctsrf(i,nsrf)
3941!
3942              zxfluxt_w(i,k) = zxfluxt_w(i,k) + flux_t_w(i,k,nsrf) * pctsrf(i,nsrf)
3943              zxfluxq_w(i,k) = zxfluxq_w(i,k) + flux_q_w(i,k,nsrf) * pctsrf(i,nsrf)
3944              zxfluxu_w(i,k) = zxfluxu_w(i,k) + flux_u_w(i,k,nsrf) * pctsrf(i,nsrf)
3945              zxfluxv_w(i,k) = zxfluxv_w(i,k) + flux_v_w(i,k,nsrf) * pctsrf(i,nsrf)
3946#ifdef ISO
3947              DO ixt=1,ntraciso
3948                zxfluxxt_x(ixt,i,k) = zxfluxxt_x(ixt,i,k) + flux_xt_x(ixt,i,k,nsrf) * pctsrf(i,nsrf)
3949                zxfluxxt_w(ixt,i,k) = zxfluxxt_w(ixt,i,k) + flux_xt_w(ixt,i,k,nsrf) * pctsrf(i,nsrf)
3950              ENDDO ! DO ixt=1,ntraciso
3951#endif
3952            ENDDO
3953          ENDDO
3954        ENDDO
3955
3956    DO i = 1, klon
3957      zxsens_x(i) = - zxfluxt_x(i,1)
3958      zxsens_w(i) = - zxfluxt_w(i,1)
3959    ENDDO
3960!!!
3961       ENDIF  ! (iflag_split .ge.1)
3962!!!
3963
3964    DO nsrf = 1, nbsrf
3965       DO k = 1, klev
3966          DO i = 1, klon
3967             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
3968             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
3969             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
3970             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
3971#ifdef ISO
3972             DO ixt=1,niso
3973               zxfluxxt(ixt,i,k) = zxfluxxt(ixt,i,k) + flux_xt(ixt,i,k,nsrf) * pctsrf(i,nsrf)
3974             ENDDO ! DO ixt=1,niso
3975#endif
3976          ENDDO
3977       ENDDO
3978    ENDDO
3979
3980    DO i = 1, klon
3981       zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
3982       zxevap(i)     = - zxfluxq(i,1) ! flux d'evaporation au sol
3983       fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
3984    ENDDO
3985
3986    ! if blowing snow
3987    if (ok_bs) then 
3988       DO nsrf = 1, nbsrf
3989       DO k = 1, klev
3990       DO i = 1, klon
3991         zxfluxqbs(i,k) = zxfluxqbs(i,k) + flux_qbs(i,k,nsrf) * pctsrf(i,nsrf)
3992       ENDDO
3993       ENDDO
3994       ENDDO
3995
3996       DO i = 1, klon
3997        zxsnowerosion(i)     = zxfluxqbs(i,1) ! blowings snow flux at the surface
3998       END DO
3999    endif
4000
4001#ifdef ISO
4002    DO i = 1, klon
4003      DO ixt=1,ntraciso
4004        zxxtevap(ixt,i)     = - zxfluxxt(ixt,i,1)
4005      ENDDO
4006    ENDDO
4007#endif
4008
4009!!!
4010
4011!
4012! Incrementer la temperature du sol
4013!
4014    zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
4015    zt2m(:) = 0.0    ; zq2m(:) = 0.0 ; zn2mout(:,:) = 0
4016    zustar(:)=0.0 ; zu10m(:) = 0.0   ; zv10m(:) = 0.0
4017    s_pblh(:) = 0.0  ; s_plcl(:) = 0.0
4018!!! jyg le 07/02/2012
4019     s_pblh_x(:) = 0.0  ; s_plcl_x(:) = 0.0
4020     s_pblh_w(:) = 0.0  ; s_plcl_w(:) = 0.0
4021!!!
4022    s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
4023    s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
4024    s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
4025    s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
4026    wstar(:,is_ave)=0.
4027   
4028!   print*,'OK pbl 9'
4029   
4030!!! nrlmd le 02/05/2011
4031    zxfluxlat_x(:) = 0.0  ;  zxfluxlat_w(:) = 0.0
4032!!!
4033   
4034    DO nsrf = 1, nbsrf
4035       DO i = 1, klon         
4036          ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
4037         
4038          wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
4039               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
4040
4041          wfevap(i,nsrf) = evap(i,nsrf)*pctsrf(i,nsrf)
4042
4043          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf(i,nsrf)
4044          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
4045       ENDDO
4046    ENDDO
4047!
4048!<al1 order 2 correction to zxtsol, for radiation computations (main atm effect of Ts)
4049   IF (iflag_order2_sollw == 1) THEN
4050    meansqT(:) = 0. ! as working buffer
4051    DO nsrf = 1, nbsrf
4052     DO i = 1, klon
4053      meansqT(i) = meansqT(i)+(ts(i,nsrf)-zxtsol(i))**2 *pctsrf(i,nsrf)
4054     ENDDO
4055    ENDDO
4056    zxtsol(:) = zxtsol(:)+1.5*meansqT(:)/zxtsol(:)
4057   ENDIF   ! iflag_order2_sollw == 1
4058!>al1
4059         
4060!!! jyg le 07/02/2012
4061       IF (iflag_split .eq.0) THEN
4062        DO nsrf = 1, nbsrf
4063         DO i = 1, klon         
4064          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
4065          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
4066!
4067          DO k = 1, 6
4068           zn2mout(i,k)  = zn2mout(i,k)  + n2mout(i,nsrf,k)  * pctsrf(i,nsrf)
4069          ENDDO 
4070!
4071          zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
4072          wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
4073          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
4074          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
4075
4076          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf(i,nsrf)
4077          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf(i,nsrf)
4078          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf(i,nsrf)
4079          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
4080          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
4081          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf(i,nsrf)
4082          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf(i,nsrf)
4083          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf(i,nsrf)
4084          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf(i,nsrf)
4085          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf(i,nsrf)
4086         ENDDO
4087        ENDDO
4088       ELSE  !(iflag_split .eq.0)
4089        DO nsrf = 1, nbsrf
4090         DO i = 1, klon         
4091!!! nrlmd le 02/05/2011
4092          zxfluxlat_x(i) = zxfluxlat_x(i) + fluxlat_x(i,nsrf) * pctsrf(i,nsrf)
4093          zxfluxlat_w(i) = zxfluxlat_w(i) + fluxlat_w(i,nsrf) * pctsrf(i,nsrf)
4094!!!
4095!!! jyg le 08/02/2012
4096!!  Pour le moment, on sort les valeurs dans (x) et (w) de pblh et de plcl ;
4097!!  pour zt2m, on fait la moyenne surfacique sur les sous-surfaces ;
4098!!  pour qsat2m, on fait la moyenne surfacique sur (x) et (w) ;
4099!!  pour les autres variables, on sort les valeurs de la region (x).
4100          zt2m(i)  = zt2m(i)  + (t2m_x(i,nsrf)+wake_s(i)*(t2m_w(i,nsrf)-t2m_x(i,nsrf))) * pctsrf(i,nsrf)
4101          zq2m(i)  = zq2m(i)  + q2m_x(i,nsrf)  * pctsrf(i,nsrf)
4102!
4103          DO k = 1, 6
4104           zn2mout(i,k)  = zn2mout(i,k)  + n2mout_x(i,nsrf,k)  * pctsrf(i,nsrf)
4105          ENDDO
4106!
4107          zustar(i) = zustar(i) + ustar_x(i,nsrf) * pctsrf(i,nsrf)
4108          wstar(i,is_ave)=wstar(i,is_ave)+wstar_x(i,nsrf)*pctsrf(i,nsrf)
4109          zu10m(i) = zu10m(i) + u10m_x(i,nsrf) * pctsrf(i,nsrf)
4110          zv10m(i) = zv10m(i) + v10m_x(i,nsrf) * pctsrf(i,nsrf)
4111!
4112          s_pblh(i)     = s_pblh(i)     + pblh_x(i,nsrf)  * pctsrf(i,nsrf)
4113          s_pblh_x(i)   = s_pblh_x(i)   + pblh_x(i,nsrf)  * pctsrf(i,nsrf)
4114          s_pblh_w(i)   = s_pblh_w(i)   + pblh_w(i,nsrf)  * pctsrf(i,nsrf)
4115!
4116          s_plcl(i)     = s_plcl(i)     + plcl_x(i,nsrf)  * pctsrf(i,nsrf)
4117          s_plcl_x(i)   = s_plcl_x(i)   + plcl_x(i,nsrf)  * pctsrf(i,nsrf)
4118          s_plcl_w(i)   = s_plcl_w(i)   + plcl_w(i,nsrf)  * pctsrf(i,nsrf)
4119!
4120          s_capCL(i)  = s_capCL(i)  + capCL_x(i,nsrf) * pctsrf(i,nsrf)
4121          s_oliqCL(i) = s_oliqCL(i) + oliqCL_x(i,nsrf)* pctsrf(i,nsrf)
4122          s_cteiCL(i) = s_cteiCL(i) + cteiCL_x(i,nsrf)* pctsrf(i,nsrf)
4123          s_pblT(i)   = s_pblT(i)   + pblT_x(i,nsrf)  * pctsrf(i,nsrf)
4124          s_therm(i)  = s_therm(i)  + therm_x(i,nsrf) * pctsrf(i,nsrf)
4125          s_trmb1(i)  = s_trmb1(i)  + trmb1_x(i,nsrf) * pctsrf(i,nsrf)
4126          s_trmb2(i)  = s_trmb2(i)  + trmb2_x(i,nsrf) * pctsrf(i,nsrf)
4127          s_trmb3(i)  = s_trmb3(i)  + trmb3_x(i,nsrf) * pctsrf(i,nsrf)
4128         ENDDO
4129        ENDDO
4130        DO i = 1, klon         
4131          qsat2m(i)= qsat2m_x(i)+ wake_s(i)*(qsat2m_x(i)-qsat2m_w(i))
4132        ENDDO
4133!!!
4134       ENDIF  ! (iflag_split .eq.0)
4135!!!
4136
4137    IF (check) THEN
4138       amn=MIN(ts(1,is_ter),1000.)
4139       amx=MAX(ts(1,is_ter),-1000.)
4140       DO i=2, klon
4141          amn=MIN(ts(i,is_ter),amn)
4142          amx=MAX(ts(i,is_ter),amx)
4143       ENDDO
4144       PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
4145    ENDIF
4146
4147!jg ?
4148!!$!
4149!!$! If a sub-surface does not exsist for a grid point, the mean value for all
4150!!$! sub-surfaces is distributed.
4151!!$!
4152!!$    DO nsrf = 1, nbsrf
4153!!$       DO i = 1, klon
4154!!$          IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
4155!!$             ts(i,nsrf)     = zxtsol(i)
4156!!$             t2m(i,nsrf)    = zt2m(i)
4157!!$             q2m(i,nsrf)    = zq2m(i)
4158!!$             u10m(i,nsrf)   = zu10m(i)
4159!!$             v10m(i,nsrf)   = zv10m(i)
4160!!$
4161!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
4162!!$             pblh(i,nsrf)   = s_pblh(i)
4163!!$             plcl(i,nsrf)   = s_plcl(i)
4164!!$             capCL(i,nsrf)  = s_capCL(i)
4165!!$             oliqCL(i,nsrf) = s_oliqCL(i)
4166!!$             cteiCL(i,nsrf) = s_cteiCL(i)
4167!!$             pblT(i,nsrf)   = s_pblT(i)
4168!!$             therm(i,nsrf)  = s_therm(i)
4169!!$             trmb1(i,nsrf)  = s_trmb1(i)
4170!!$             trmb2(i,nsrf)  = s_trmb2(i)
4171!!$             trmb3(i,nsrf)  = s_trmb3(i)
4172!!$          ENDIF
4173!!$       ENDDO
4174!!$    ENDDO
4175
4176
4177    DO i = 1, klon
4178       fder(i) = - 4.0*RSIGMA*zxtsol(i)**3
4179    ENDDO
4180   
4181    zxqsurf(:) = 0.0
4182    zxsnow(:)  = 0.0
4183#ifdef ISO
4184    zxxtsnow(:,:)  = 0.0
4185#endif
4186
4187    DO nsrf = 1, nbsrf
4188       DO i = 1, klon
4189          zxqsurf(i) = zxqsurf(i) + MAX(qsurf(i,nsrf),0.0) * pctsrf(i,nsrf)
4190          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
4191#ifdef ISO
4192          DO ixt=1,niso
4193            zxxtsnow(ixt,i)  = zxxtsnow(ixt,i)  + xtsnow(ixt,i,nsrf)  * pctsrf(i,nsrf)
4194          ENDDO ! DO ixt=1,niso
4195#endif
4196       ENDDO
4197    ENDDO
4198
4199! Premier niveau de vent sortie dans physiq.F
4200    zu1(:) = u(:,1)
4201    zv1(:) = v(:,1)
4202
4203  END SUBROUTINE pbl_surface
4204!
4205!****************************************************************************************
4206!
4207  SUBROUTINE pbl_surface_final(fder_rst, snow_rst, qsurf_rst, ftsoil_rst &
4208#ifdef ISO
4209       ,xtsnow_rst,Rland_ice_rst &
4210#endif       
4211       )
4212
4213    USE indice_sol_mod
4214#ifdef ISO
4215#ifdef ISOVERIF
4216    USE isotopes_mod, ONLY: iso_eau,ridicule
4217    USE isotopes_verif_mod, ONLY: errmax,errmaxrel
4218#endif   
4219#endif
4220    USE dimsoil_mod_h, ONLY: nsoilmx
4221
4222! Ouput variables
4223!****************************************************************************************
4224    REAL, DIMENSION(klon), INTENT(OUT)                 :: fder_rst
4225    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: snow_rst
4226    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
4227    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
4228#ifdef ISO
4229    REAL, DIMENSION(niso,klon, nbsrf), INTENT(OUT)     :: xtsnow_rst
4230    REAL, DIMENSION(niso,klon), INTENT(OUT)            :: Rland_ice_rst
4231#endif
4232
4233 
4234!****************************************************************************************
4235! Return module variables for writing to restart file
4236!
4237!****************************************************************************************   
4238    fder_rst(:)       = fder(:)
4239    snow_rst(:,:)     = snow(:,:)
4240    qsurf_rst(:,:)    = qsurf(:,:)
4241    ftsoil_rst(:,:,:) = ftsoil(:,:,:)
4242#ifdef ISO
4243    xtsnow_rst(:,:,:)  = xtsnow(:,:,:)
4244    Rland_ice_rst(:,:) = Rland_ice(:,:)
4245#endif
4246
4247!****************************************************************************************
4248! Deallocate module variables
4249!
4250!****************************************************************************************
4251!   DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
4252    IF (ALLOCATED(fder)) DEALLOCATE(fder)
4253    IF (ALLOCATED(snow)) DEALLOCATE(snow)
4254    IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
4255    IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
4256    IF (ALLOCATED(ydTs0)) DEALLOCATE(ydTs0)
4257    IF (ALLOCATED(ydqs0)) DEALLOCATE(ydqs0)
4258#ifdef ISO
4259    IF (ALLOCATED(xtsnow)) DEALLOCATE(xtsnow)
4260    IF (ALLOCATED(Rland_ice)) DEALLOCATE(Rland_ice)
4261    IF (ALLOCATED(Roce)) DEALLOCATE(Roce)
4262#endif
4263
4264!jyg<
4265!****************************************************************************************
4266! Deallocate variables for pbl splitting
4267!
4268!****************************************************************************************
4269
4270    CALL wx_pbl_final
4271!>jyg
4272
4273  END SUBROUTINE pbl_surface_final
4274
4275!****************************************************************************************
4276!
4277
4278!albedo SB >>>
4279  SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, &
4280       evap, z0m, z0h, agesno,                                  &
4281       tsurf,alb_dir,alb_dif, ustar, u10m, v10m, tke &
4282#ifdef ISO
4283      ,xtevap  &
4284#endif
4285&      ) 
4286    !albedo SB <<<
4287    ! Give default values where new fraction has appread
4288
4289USE compbl_mod_h
4290        USE clesphys_mod_h
4291    USE indice_sol_mod
4292    use phys_state_var_mod, only: delta_sal, ds_ns, dt_ns, delta_sst, dter, &
4293         dser, dt_ds
4294    use config_ocean_skin_m, only: activate_ocean_skin
4295
4296! Input variables
4297!****************************************************************************************
4298    INTEGER, INTENT(IN)                     :: itime
4299    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
4300
4301! InOutput variables
4302!****************************************************************************************
4303    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: tsurf
4304!albedo SB >>>
4305    REAL, DIMENSION(klon,nsw,nbsrf), INTENT(INOUT)       :: alb_dir, alb_dif
4306    INTEGER :: k
4307!albedo SB <<<
4308    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: ustar,u10m, v10m
4309    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: evap, agesno
4310    REAL, DIMENSION(klon,nbsrf+1), INTENT(INOUT)        :: z0m,z0h
4311    REAL, DIMENSION(klon,klev+1,nbsrf+1), INTENT(INOUT) :: tke
4312#ifdef ISO
4313    REAL, DIMENSION(ntraciso,klon,nbsrf), INTENT(INOUT)        :: xtevap
4314#endif
4315
4316! Local variables
4317!****************************************************************************************
4318    INTEGER           :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
4319    CHARACTER(len=80) :: abort_message
4320    CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
4321    INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
4322#ifdef ISO
4323    INTEGER           :: ixt
4324#endif
4325!
4326! All at once !!
4327!****************************************************************************************
4328   
4329    DO nsrf = 1, nbsrf
4330       ! First decide complement sub-surfaces
4331       SELECT CASE (nsrf)
4332       CASE(is_oce)
4333          nsrf_comp1=is_sic
4334          nsrf_comp2=is_ter
4335          nsrf_comp3=is_lic
4336       CASE(is_sic)
4337          nsrf_comp1=is_oce
4338          nsrf_comp2=is_ter
4339          nsrf_comp3=is_lic
4340       CASE(is_ter)
4341          nsrf_comp1=is_lic
4342          nsrf_comp2=is_oce
4343          nsrf_comp3=is_sic
4344       CASE(is_lic)
4345          nsrf_comp1=is_ter
4346          nsrf_comp2=is_oce
4347          nsrf_comp3=is_sic
4348       END SELECT
4349
4350       ! Initialize all new fractions
4351       DO i=1, klon
4352          IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
4353             
4354             IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
4355                ! Use the complement sub-surface, keeping the continents unchanged
4356                qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
4357                evap(i,nsrf)  = evap(i,nsrf_comp1)
4358                z0m(i,nsrf) = z0m(i,nsrf_comp1)
4359                z0h(i,nsrf) = z0h(i,nsrf_comp1)
4360                tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
4361!albedo SB >>>
4362                DO k=1,nsw
4363                 alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp1)
4364                 alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp1)
4365                ENDDO
4366!albedo SB <<<
4367                ustar(i,nsrf)  = ustar(i,nsrf_comp1)
4368                u10m(i,nsrf)  = u10m(i,nsrf_comp1)
4369                v10m(i,nsrf)  = v10m(i,nsrf_comp1)
4370#ifdef ISO
4371                DO ixt=1,ntraciso
4372                  xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp1)       
4373                ENDDO       
4374#endif
4375                IF (iflag_pbl > 1) THEN
4376                 tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
4377                ENDIF
4378                mfois(nsrf) = mfois(nsrf) + 1
4379                ! F. Codron sensible default values for ocean and sea ice
4380                IF (nsrf.EQ.is_oce) THEN
4381                   tsurf(i,nsrf) = 271.35
4382                   ! (temperature of sea water under sea ice, so that
4383                   ! is also the temperature of appearing sea water)
4384                   DO k=1,nsw
4385                      alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
4386                      alb_dif(i,k,nsrf) = 0.06
4387                   ENDDO
4388                   if (activate_ocean_skin >= 1) then
4389                      if (activate_ocean_skin == 2 &
4390                           .and. type_ocean == "couple") then
4391                         delta_sal(i) = 0.
4392                         delta_sst(i) = 0.
4393                         dter(i) = 0.
4394                         dser(i) = 0.
4395                         dt_ds(i) = 0.
4396                      end if
4397                     
4398                      ds_ns(i) = 0.
4399                      dt_ns(i) = 0.
4400                   end if
4401                ELSE IF (nsrf.EQ.is_sic) THEN
4402                   tsurf(i,nsrf) = 271.35
4403                   ! (Temperature at base of sea ice. Surface
4404                   ! temperature could be higher, up to 0 Celsius
4405                   ! degrees. We set it to -1.8 Celsius degrees for
4406                   ! consistency with the ocean slab model.)
4407                   DO k=1,nsw
4408                      alb_dir(i,k,nsrf) = 0.3 ! thin ice
4409                      alb_dif(i,k,nsrf) = 0.3
4410                   ENDDO
4411                ENDIF
4412             ELSE
4413                ! The continents have changed. The new fraction receives the mean sum of the existent fractions
4414                qsurf(i,nsrf) = qsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
4415                evap(i,nsrf)  = evap(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
4416                z0m(i,nsrf) = z0m(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + z0m(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
4417                z0h(i,nsrf) = z0h(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + z0h(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
4418                tsurf(i,nsrf) = tsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
4419!albedo SB >>>
4420                DO k=1,nsw
4421                 alb_dir(i,k,nsrf)=alb_dir(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
4422                                        alb_dir(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
4423                 alb_dif(i,k,nsrf)=alb_dif(i,k,nsrf_comp2)*pctsrf_old(i,nsrf_comp2)+&
4424                                        alb_dif(i,k,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
4425                ENDDO
4426!albedo SB <<<
4427                ustar(i,nsrf)  = ustar(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
4428                u10m(i,nsrf)  = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
4429                v10m(i,nsrf)  = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
4430#ifdef ISO
4431                DO ixt=1,ntraciso
4432                  xtevap(ixt,i,nsrf) = xtevap(ixt,i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) &
4433                                     + xtevap(ixt,i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
4434                ENDDO       
4435#endif
4436                IF (iflag_pbl > 1) THEN
4437                 tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
4438                ENDIF
4439           
4440                ! Security abort. This option has never been tested. To test, comment the following line.
4441!                abort_message='The fraction of the continents have changed!'
4442!                CALL abort_physic(modname,abort_message,1)
4443                nfois(nsrf) = nfois(nsrf) + 1
4444             ENDIF
4445             snow(i,nsrf)     = 0.
4446             agesno(i,nsrf)   = 0.
4447             ftsoil(i,:,nsrf) = tsurf(i,nsrf)
4448#ifdef ISO           
4449             xtsnow(:,i,nsrf) = 0.
4450#endif
4451          ELSE
4452             pfois(nsrf) = pfois(nsrf)+ 1
4453          ENDIF
4454       ENDDO
4455       
4456    ENDDO
4457
4458  END SUBROUTINE pbl_surface_newfrac
4459
4460!****************************************************************************************
4461
4462END MODULE pbl_surface_mod
Note: See TracBrowser for help on using the repository browser.