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

Last change on this file since 5314 was 2952, checked in by Laurent Fairhead, 7 years ago

Parametrization of drag by copses
Need version 4465 of ORCHIDEE at least

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