source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phylmd/pbl_surface_mod.F90 @ 3603

Last change on this file since 3603 was 3406, checked in by jghattas, 6 years ago

Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore

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