source: LMDZ6/branches/LMDZ-tracers/libf/phylmd/pbl_surface_mod.F90 @ 3957

Last change on this file since 3957 was 3851, checked in by dcugnet, 4 years ago

Update the branch to the current trunk.

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