source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/pbl_surface_mod.F90 @ 5403

Last change on this file since 5403 was 4727, checked in by idelkadi, 14 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch

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