source: LMDZ5/branches/IPSLCM5A2.1/libf/phylmd/pbl_surface_mod.F90 @ 5017

Last change on this file since 5017 was 2455, checked in by jyg, 9 years ago

Implementation of a second order distribution on

sub-surfaces of longwave net radiance (Alain
Lahellec).

Introduction of the flag iflag_order2_sollw:

=1 ==> second order distribution.
=0 ==> linear distribution (Default).

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