source: LMDZ5/trunk/libf/phylmd/pbl_surface_mod.F90 @ 2345

Last change on this file since 2345 was 2345, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • move test_disvert_m to dynlonlat_phylonlat/phylmd since it is only used by ce0l and relies on dynamics.
  • put "config_inca" in tracinca_mod so physics routines can get the info from there rather than from control_mod.
  • get rid of references to "control_mod" from within the physics.

EM

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