My Project
 All Classes Files Functions Variables Macros
pbl_surface_mod.F90
Go to the documentation of this file.
1 !
2 ! $Id: pbl_surface_mod.F90 1761 2013-06-04 10:11:38Z fairhead $
3 !
4 MODULE pbl_surface_mod
5 !
6 ! Planetary Boundary Layer and Surface module
7 !
8 ! This module manage the calculation of turbulent diffusion in the boundary layer
9 ! and all interactions towards the differents sub-surfaces.
10 !
11 !
12  USE dimphy
13  USE mod_phys_lmdz_para, ONLY : mpi_size
14  USE ioipsl
15  USE surface_data, ONLY : type_ocean, ok_veget
16  USE surf_land_mod, ONLY : surf_land
17  USE surf_landice_mod, ONLY : surf_landice
18  USE surf_ocean_mod, ONLY : surf_ocean
19  USE surf_seaice_mod, ONLY : surf_seaice
20  USE cpl_mod, ONLY : gath2cpl
21  USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up
22  USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up
23  USE coef_diff_turb_mod, ONLY : coef_diff_turb
24  USE control_mod
25 
26 
27  IMPLICIT NONE
28 
29 ! Declaration of variables saved in restart file
30  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: qsol ! water height in the soil (mm)
31  !$OMP THREADPRIVATE(qsol)
32  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: fder ! flux drift
33  !$OMP THREADPRIVATE(fder)
34  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: snow ! snow at surface
35  !$OMP THREADPRIVATE(snow)
36  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qsurf ! humidity at surface
37  !$OMP THREADPRIVATE(qsurf)
38  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: evap ! evaporation at surface
39  !$OMP THREADPRIVATE(evap)
40  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: rugos ! rugosity at surface (m)
41  !$OMP THREADPRIVATE(rugos)
42  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: agesno ! age of snow at surface
43  !$OMP THREADPRIVATE(agesno)
44  REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: ftsoil ! soil temperature
45  !$OMP THREADPRIVATE(ftsoil)
46 
47 CONTAINS
48 !
49 !****************************************************************************************
50 !
51  SUBROUTINE pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst,&
52  evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
53 
54 ! This routine should be called after the restart file has been read.
55 ! This routine initialize the restart variables and does some validation tests
56 ! for the index of the different surfaces and tests the choice of type of ocean.
57 
58  include "indicesol.h"
59  include "dimsoil.h"
60  include "iniprint.h"
61 
62 ! Input variables
63 !****************************************************************************************
64  REAL, DIMENSION(klon), INTENT(IN) :: qsol_rst
65  REAL, DIMENSION(klon), INTENT(IN) :: fder_rst
66  REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: snow_rst
67  REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: qsurf_rst
68  REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: evap_rst
69  REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: rugos_rst
70  REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: agesno_rst
71  REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
72 
73 
74 ! Local variables
75 !****************************************************************************************
76  INTEGER :: ierr
77  CHARACTER(len=80) :: abort_message
78  CHARACTER(len = 20) :: modname = 'pbl_surface_init'
79 
80 
81 !****************************************************************************************
82 ! Allocate and initialize module variables with fields read from restart file.
83 !
84 !****************************************************************************************
85  ALLOCATE(qsol(klon), stat=ierr)
86  IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
87 
88  ALLOCATE(fder(klon), stat=ierr)
89  IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
90 
91  ALLOCATE(snow(klon,nbsrf), stat=ierr)
92  IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
93 
94  ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
95  IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
96 
97  ALLOCATE(evap(klon,nbsrf), stat=ierr)
98  IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
99 
100  ALLOCATE(rugos(klon,nbsrf), stat=ierr)
101  IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
102 
103  ALLOCATE(agesno(klon,nbsrf), stat=ierr)
104  IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
105 
106  ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
107  IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
108 
109 
110  qsol(:) = qsol_rst(:)
111  fder(:) = fder_rst(:)
112  snow(:,:) = snow_rst(:,:)
113  qsurf(:,:) = qsurf_rst(:,:)
114  evap(:,:) = evap_rst(:,:)
115  rugos(:,:) = rugos_rst(:,:)
116  agesno(:,:) = agesno_rst(:,:)
117  ftsoil(:,:,:) = ftsoil_rst(:,:,:)
118 
119 
120 !****************************************************************************************
121 ! Test for sub-surface indices
122 !
123 !****************************************************************************************
124  IF (is_ter /= 1) THEN
125  WRITE(lunout,*)" *** Warning ***"
126  WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
127  WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
128  abort_message="voir ci-dessus"
129  CALL abort_gcm(modname,abort_message,1)
130  ENDIF
131 
132  IF ( is_oce > is_sic ) THEN
133  WRITE(lunout,*)' *** Warning ***'
134  WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
135  WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
136  WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
137  abort_message='voir ci-dessus'
138  CALL abort_gcm(modname,abort_message,1)
139  ENDIF
140 
141  IF ( is_lic > is_sic ) THEN
142  WRITE(lunout,*)' *** Warning ***'
143  WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
144  WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
145  WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
146  abort_message='voir ci-dessus'
147  CALL abort_gcm(modname,abort_message,1)
148  ENDIF
149 
150 !****************************************************************************************
151 ! Validation of ocean mode
152 !
153 !****************************************************************************************
154 
155  IF (type_ocean /= 'slab ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
156  WRITE(lunout,*)' *** Warning ***'
157  WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
158  abort_message='option pour l''ocean non valable'
159  CALL abort_gcm(modname,abort_message,1)
160  ENDIF
161 
162  END SUBROUTINE pbl_surface_init
163 !
164 !****************************************************************************************
165 !
166 
167  SUBROUTINE pbl_surface( &
168  dtime, date0, itap, jour, &
169  debut, lafin, &
170  rlon, rlat, rugoro, rmu0, &
171  rain_f, snow_f, solsw_m, sollw_m, &
172  t, q, u, v, &
173  pplay, paprs, pctsrf, &
174  ts, alb1, alb2,ustar, u10m, v10m, &
175  lwdown_m, cdragh, cdragm, zu1, zv1, &
176  alb1_m, alb2_m, zxsens, zxevap, &
177  zxtsol, zxfluxlat, zt2m, qsat2m, &
178  d_t, d_q, d_u, d_v, d_t_diss, &
179  zcoefh, zcoefm, slab_wfbils, &
180  qsol_d, zq2m, s_pblh, s_plcl, &
181  s_capcl, s_oliqcl, s_cteicl, s_pblt, &
182  s_therm, s_trmb1, s_trmb2, s_trmb3, &
183  zxrugs,zustar,zu10m, zv10m, fder_print, &
184  zxqsurf, rh2m, zxfluxu, zxfluxv, &
185  rugos_d, agesno_d, sollw, solsw, &
186  d_ts, evap_d, fluxlat, t2m, &
187  wfbils, wfbilo, flux_t, flux_u, flux_v,&
188  dflux_t, dflux_q, zxsnow, &
189  zxfluxt, zxfluxq, q2m, flux_q, tke )
190 !****************************************************************************************
191 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
192 ! Objet: interface de "couche limite" (diffusion verticale)
193 !
194 !AA REM:
195 !AA-----
196 !AA Tout ce qui a trait au traceurs est dans phytrac maintenant
197 !AA pour l'instant le calcul de la couche limite pour les traceurs
198 !AA se fait avec cltrac et ne tient pas compte de la differentiation
199 !AA des sous-fraction de sol.
200 !AA REM bis :
201 !AA----------
202 !AA Pour pouvoir extraire les coefficient d'echanges et le vent
203 !AA dans la premiere couche, 3 champs supplementaires ont ete crees
204 !AA zcoefh, zu1 et zv1. Pour l'instant nous avons moyenne les valeurs
205 !AA de ces trois champs sur les 4 subsurfaces du modele. Dans l'avenir
206 !AA si les informations des subsurfaces doivent etre prises en compte
207 !AA il faudra sortir ces memes champs en leur ajoutant une dimension,
208 !AA c'est a dire nbsrf (nbre de subsurface).
209 !
210 ! Arguments:
211 !
212 ! dtime----input-R- interval du temps (secondes)
213 ! itap-----input-I- numero du pas de temps
214 ! date0----input-R- jour initial
215 ! t--------input-R- temperature (K)
216 ! q--------input-R- vapeur d'eau (kg/kg)
217 ! u--------input-R- vitesse u
218 ! v--------input-R- vitesse v
219 ! ts-------input-R- temperature du sol (en Kelvin)
220 ! paprs----input-R- pression a intercouche (Pa)
221 ! pplay----input-R- pression au milieu de couche (Pa)
222 ! rlat-----input-R- latitude en degree
223 ! rugos----input-R- longeur de rugosite (en m)
224 !
225 ! d_t------output-R- le changement pour "t"
226 ! d_q------output-R- le changement pour "q"
227 ! d_u------output-R- le changement pour "u"
228 ! d_v------output-R- le changement pour "v"
229 ! d_ts-----output-R- le changement pour "ts"
230 ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
231 ! (orientation positive vers le bas)
232 ! tke---input/output-R- tke (kg/m**2/s)
233 ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
234 ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
235 ! flux_v---output-R- tension du vent Y: (kg m/s)/(m**2 s) ou Pascal
236 ! dflux_t--output-R- derive du flux sensible
237 ! dflux_q--output-R- derive du flux latent
238 ! zu1------output-R- le vent dans la premiere couche
239 ! zv1------output-R- le vent dans la premiere couche
240 ! trmb1----output-R- deep_cape
241 ! trmb2----output-R- inhibition
242 ! trmb3----output-R- Point Omega
243 ! cteiCL---output-R- Critere d'instab d'entrainmt des nuages de CL
244 ! plcl-----output-R- Niveau de condensation
245 ! pblh-----output-R- HCL
246 ! pblT-----output-R- T au nveau HCL
247 !
248  USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
249  IMPLICIT NONE
250 
251  include "indicesol.h"
252  include "dimsoil.h"
253  include "YOMCST.h"
254  include "iniprint.h"
255  include "FCTTRE.h"
256  include "clesphys.h"
257  include "compbl.h"
258  include "dimensions.h"
259  include "YOETHF.h"
260  include "temps.h"
261 !****************************************************************************************
262 ! Declarations specifiques pour le 1D. A reprendre
263 ! Input variables
264 !****************************************************************************************
265  REAL, INTENT(IN) :: dtime ! time interval (s)
266  REAL, INTENT(IN) :: date0 ! initial day
267  INTEGER, INTENT(IN) :: itap ! time step
268  INTEGER, INTENT(IN) :: jour ! current day of the year
269  LOGICAL, INTENT(IN) :: debut ! true if first run step
270  LOGICAL, INTENT(IN) :: lafin ! true if last run step
271  REAL, DIMENSION(klon), INTENT(IN) :: rlon ! longitudes in degrees
272  REAL, DIMENSION(klon), INTENT(IN) :: rlat ! latitudes in degrees
273  REAL, DIMENSION(klon), INTENT(IN) :: rugoro ! rugosity length
274  REAL, DIMENSION(klon), INTENT(IN) :: rmu0 ! cosine of solar zenith angle
275  REAL, DIMENSION(klon), INTENT(IN) :: rain_f ! rain fall
276  REAL, DIMENSION(klon), INTENT(IN) :: snow_f ! snow fall
277  REAL, DIMENSION(klon), INTENT(IN) :: solsw_m ! net shortwave radiation at mean surface
278  REAL, DIMENSION(klon), INTENT(IN) :: sollw_m ! net longwave radiation at mean surface
279  REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K)
280  REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! water vapour (kg/kg)
281  REAL, DIMENSION(klon,klev), INTENT(IN) :: u ! u speed
282  REAL, DIMENSION(klon,klev), INTENT(IN) :: v ! v speed
283  REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pression (Pa)
284  REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression between layers (Pa)
285  REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! sub-surface fraction
286 
287 ! Input/Output variables
288 !****************************************************************************************
289  REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ts ! temperature at surface (K)
290  REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: alb1 ! albedo in visible SW interval
291  REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: alb2 ! albedo in near infra-red SW interval
292  REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: ustar ! u* (m/s)
293  REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: u10m ! u speed at 10m
294  REAL, DIMENSION(klon, nbsrf), INTENT(INOUT) :: v10m ! v speed at 10m
295  REAL, DIMENSION(klon, klev+1, nbsrf+1), INTENT(INOUT) :: tke
296 
297 ! Output variables
298 !****************************************************************************************
299  REAL, DIMENSION(klon), INTENT(OUT) :: lwdown_m ! Downcoming longwave radiation
300  REAL, DIMENSION(klon), INTENT(OUT) :: cdragh ! drag coefficient for T and Q
301  REAL, DIMENSION(klon), INTENT(OUT) :: cdragm ! drag coefficient for wind
302  REAL, DIMENSION(klon), INTENT(OUT) :: zu1 ! u wind speed in first layer
303  REAL, DIMENSION(klon), INTENT(OUT) :: zv1 ! v wind speed in first layer
304  REAL, DIMENSION(klon), INTENT(OUT) :: alb1_m ! mean albedo in visible SW interval
305  REAL, DIMENSION(klon), INTENT(OUT) :: alb2_m ! mean albedo in near IR SW interval
306  REAL, DIMENSION(klon), INTENT(OUT) :: zxsens ! sensible heat flux at surface with inversed sign
307  ! (=> positive sign upwards)
308  REAL, DIMENSION(klon), INTENT(OUT) :: zxevap ! water vapour flux at surface, positiv upwards
309  REAL, DIMENSION(klon), INTENT(OUT) :: zxtsol ! temperature at surface, mean for each grid point
310  REAL, DIMENSION(klon), INTENT(OUT) :: zxfluxlat ! latent flux, mean for each grid point
311  REAL, DIMENSION(klon), INTENT(OUT) :: zt2m ! temperature at 2m, mean for each grid point
312  REAL, DIMENSION(klon), INTENT(OUT) :: qsat2m
313  REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t ! change in temperature
314  REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_diss ! change in temperature
315  REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_q ! change in water vapour
316  REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u ! change in u speed
317  REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_v ! change in v speed
318  REAL, DIMENSION(klon, klev,nbsrf+1), INTENT(OUT) :: zcoefh ! coef for turbulent diffusion of T and Q, mean for each grid point
319  REAL, DIMENSION(klon, klev,nbsrf+1), INTENT(OUT) :: zcoefm ! coef for turbulent diffusion of U and V (?), mean for each grid point
320 
321 ! Output only for diagnostics
322  REAL, DIMENSION(klon), INTENT(OUT) :: slab_wfbils! heat balance at surface only for slab at ocean points
323  REAL, DIMENSION(klon), INTENT(OUT) :: qsol_d ! water height in the soil (mm)
324  REAL, DIMENSION(klon), INTENT(OUT) :: zq2m ! water vapour at 2m, mean for each grid point
325  REAL, DIMENSION(klon), INTENT(OUT) :: s_pblh ! height of the planetary boundary layer(HPBL)
326  REAL, DIMENSION(klon), INTENT(OUT) :: s_plcl ! condensation level
327  REAL, DIMENSION(klon), INTENT(OUT) :: s_capcl ! CAPE of PBL
328  REAL, DIMENSION(klon), INTENT(OUT) :: s_oliqcl ! liquid water intergral of PBL
329  REAL, DIMENSION(klon), INTENT(OUT) :: s_cteicl ! cloud top instab. crit. of PBL
330  REAL, DIMENSION(klon), INTENT(OUT) :: s_pblt ! temperature at PBLH
331  REAL, DIMENSION(klon), INTENT(OUT) :: s_therm ! thermal virtual temperature excess
332  REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb1 ! deep cape, mean for each grid point
333  REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb2 ! inhibition, mean for each grid point
334  REAL, DIMENSION(klon), INTENT(OUT) :: s_trmb3 ! point Omega, mean for each grid point
335  REAL, DIMENSION(klon), INTENT(OUT) :: zxrugs ! rugosity at surface (m), mean for each grid point
336  REAL, DIMENSION(klon), INTENT(OUT) :: zustar ! u*
337  REAL, DIMENSION(klon), INTENT(OUT) :: zu10m ! u speed at 10m, mean for each grid point
338  REAL, DIMENSION(klon), INTENT(OUT) :: zv10m ! v speed at 10m, mean for each grid point
339  REAL, DIMENSION(klon), INTENT(OUT) :: fder_print ! fder for printing (=fder(i) + dflux_t(i) + dflux_q(i))
340  REAL, DIMENSION(klon), INTENT(OUT) :: zxqsurf ! humidity at surface, mean for each grid point
341  REAL, DIMENSION(klon), INTENT(OUT) :: rh2m ! relative humidity at 2m
342  REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxu ! u wind tension, mean for each grid point
343  REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxv ! v wind tension, mean for each grid point
344  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_d ! rugosity length (m)
345  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_d ! age of snow at surface
346  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: solsw ! net shortwave radiation at surface
347  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: sollw ! net longwave radiation at surface
348  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: d_ts ! change in temperature at surface
349  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_d ! evaporation at surface
350  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: fluxlat ! latent flux
351  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: t2m ! temperature at 2 meter height
352  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbils ! heat balance at surface
353  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: wfbilo ! water balance at surface
354  REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_t ! sensible heat flux (CpT) J/m**2/s (W/m**2)
355  ! positve orientation downwards
356  REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_u ! u wind tension (kg m/s)/(m**2 s) or Pascal
357  REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_v ! v wind tension (kg m/s)/(m**2 s) or Pascal
358 
359 ! Output not needed
360  REAL, DIMENSION(klon), INTENT(OUT) :: dflux_t ! change of sensible heat flux
361  REAL, DIMENSION(klon), INTENT(OUT) :: dflux_q ! change of water vapour flux
362  REAL, DIMENSION(klon), INTENT(OUT) :: zxsnow ! snow at surface, mean for each grid point
363  REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxt ! sensible heat flux, mean for each grid point
364  REAL, DIMENSION(klon, klev), INTENT(OUT) :: zxfluxq ! water vapour flux, mean for each grid point
365  REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m ! water vapour at 2 meter height
366  REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q ! water vapour flux(latent flux) (kg/m**2/s)
367 
368 
369 ! Local variables with attribute SAVE
370 !****************************************************************************************
371  INTEGER, SAVE :: nhoridbg, nidbg ! variables for IOIPSL
372 !$OMP THREADPRIVATE(nhoridbg, nidbg)
373  LOGICAL, SAVE :: debugindex=.false.
374 !$OMP THREADPRIVATE(debugindex)
375  LOGICAL, SAVE :: first_call=.true.
376 !$OMP THREADPRIVATE(first_call)
377  CHARACTER(len=8), DIMENSION(nbsrf), SAVE :: cl_surf
378 !$OMP THREADPRIVATE(cl_surf)
379 
380 ! Other local variables
381 !****************************************************************************************
382  INTEGER :: i, k, nsrf
383  INTEGER :: knon, j
384  INTEGER :: idayref
385  INTEGER , DIMENSION(klon) :: ni
386  REAL :: zx_alf1, zx_alf2 !valeur ambiante par extrapola
387  REAL :: amn, amx
388  REAL :: f1 ! fraction de longeurs visibles parmi tout SW intervalle
389  REAL, DIMENSION(klon) :: r_co2_ppm ! taux CO2 atmosphere
390  REAL, DIMENSION(klon) :: yts, yrugos, ypct, yz0_new
391  REAL, DIMENSION(klon) :: yalb, yalb1, yalb2
392  REAL, DIMENSION(klon) :: yu1, yv1,ytoto
393  REAL, DIMENSION(klon) :: ysnow, yqsurf, yagesno, yqsol
394  REAL, DIMENSION(klon) :: yrain_f, ysnow_f
395  REAL, DIMENSION(klon) :: ysolsw, ysollw
396  REAL, DIMENSION(klon) :: yfder
397  REAL, DIMENSION(klon) :: yrugoro
398  REAL, DIMENSION(klon) :: yfluxlat
399  REAL, DIMENSION(klon) :: y_d_ts
400  REAL, DIMENSION(klon) :: y_flux_t1, y_flux_q1
401  REAL, DIMENSION(klon) :: y_dflux_t, y_dflux_q
402  REAL, DIMENSION(klon) :: y_flux_u1, y_flux_v1
403  REAL, DIMENSION(klon) :: yt2m, yq2m, yu10m
404  REAL, DIMENSION(klon) :: yustar
405  REAL, DIMENSION(klon) :: ywindsp
406  REAL, DIMENSION(klon) :: yt10m, yq10m
407  REAL, DIMENSION(klon) :: ypblh
408  REAL, DIMENSION(klon) :: ylcl
409  REAL, DIMENSION(klon) :: ycapcl
410  REAL, DIMENSION(klon) :: yoliqcl
411  REAL, DIMENSION(klon) :: ycteicl
412  REAL, DIMENSION(klon) :: ypblt
413  REAL, DIMENSION(klon) :: ytherm
414  REAL, DIMENSION(klon) :: ytrmb1
415  REAL, DIMENSION(klon) :: ytrmb2
416  REAL, DIMENSION(klon) :: ytrmb3
417  REAL, DIMENSION(klon) :: uzon, vmer
418  REAL, DIMENSION(klon) :: tair1, qair1, tairsol
419  REAL, DIMENSION(klon) :: psfce, patm
420  REAL, DIMENSION(klon) :: qairsol, zgeo1
421  REAL, DIMENSION(klon) :: rugo1
422  REAL, DIMENSION(klon) :: yfluxsens
423  REAL, DIMENSION(klon) :: acoefh, acoefq, bcoefh, bcoefq
424  REAL, DIMENSION(klon) :: acoefu, acoefv, bcoefu, bcoefv
425  REAL, DIMENSION(klon) :: ypsref
426  REAL, DIMENSION(klon) :: yevap, ytsurf_new, yalb1_new, yalb2_new
427  REAL, DIMENSION(klon) :: ztsol
428  REAL, DIMENSION(klon) :: alb_m ! mean albedo for whole SW interval
429  REAL, DIMENSION(klon,klev) :: y_d_t, y_d_q, y_d_t_diss
430  REAL, DIMENSION(klon,klev) :: y_d_u, y_d_v
431  REAL, DIMENSION(klon,klev) :: y_flux_t, y_flux_q
432  REAL, DIMENSION(klon,klev) :: y_flux_u, y_flux_v
433  REAL, DIMENSION(klon,klev) :: ycoefh, ycoefm,ycoefq
434  REAL, DIMENSION(klon) :: ycdragh, ycdragm
435  REAL, DIMENSION(klon,klev) :: yu, yv
436  REAL, DIMENSION(klon,klev) :: yt, yq
437  REAL, DIMENSION(klon,klev) :: ypplay, ydelp
438  REAL, DIMENSION(klon,klev) :: delp
439  REAL, DIMENSION(klon,klev+1) :: ypaprs
440  REAL, DIMENSION(klon,klev+1) :: ytke
441  REAL, DIMENSION(klon,nsoilmx) :: ytsoil
442  CHARACTER(len=80) :: abort_message
443  CHARACTER(len=20) :: modname = 'pbl_surface'
444  LOGICAL, PARAMETER :: zxli=.false. ! utiliser un jeu de fonctions simples
445  LOGICAL, PARAMETER :: check=.false.
446  REAL, DIMENSION(klon) :: kech_h ! Coefficient d'echange pour l'energie
447 
448 ! For debugging with IOIPSL
449  INTEGER, DIMENSION(iim*(jjm+1)) :: ndexbg
450  REAL :: zjulian
451  REAL, DIMENSION(klon) :: tabindx
452  REAL, DIMENSION(iim,jjm+1) :: zx_lon, zx_lat
453  REAL, DIMENSION(iim,jjm+1) :: debugtab
454 
455 
456  REAL, DIMENSION(klon,nbsrf) :: pblh ! height of the planetary boundary layer
457  REAL, DIMENSION(klon,nbsrf) :: plcl ! condensation level
458  REAL, DIMENSION(klon,nbsrf) :: capcl
459  REAL, DIMENSION(klon,nbsrf) :: oliqcl
460  REAL, DIMENSION(klon,nbsrf) :: cteicl
461  REAL, DIMENSION(klon,nbsrf) :: pblt
462  REAL, DIMENSION(klon,nbsrf) :: therm
463  REAL, DIMENSION(klon,nbsrf) :: trmb1 ! deep cape
464  REAL, DIMENSION(klon,nbsrf) :: trmb2 ! inhibition
465  REAL, DIMENSION(klon,nbsrf) :: trmb3 ! point Omega
466  REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m
467  REAL, DIMENSION(klon,nbsrf) :: zx_t1
468  REAL, DIMENSION(klon, nbsrf) :: alb ! mean albedo for whole SW interval
469  REAL, DIMENSION(klon) :: ylwdown ! jg : temporary (ysollwdown)
470 
471  REAL :: zx_qs1, zcor1, zdelta1
472 
473 !****************************************************************************************
474 ! Declarations specifiques pour le 1D. A reprendre
475 !****************************************************************************************
476  REAL :: fsens,flat
477  LOGICAL :: ok_flux_surf ! initialized during first_call below
478  COMMON /flux_arp/fsens,flat,ok_flux_surf
479 ! End of declarations
480 !****************************************************************************************
481 
482 
483 !****************************************************************************************
484 ! 1) Initialisation and validation tests
485 ! Only done first time entering this subroutine
486 !
487 !****************************************************************************************
488 
489  IF (first_call) THEN
490  first_call=.false.
491 
492  ! Initialize ok_flux_surf (for 1D model)
493  if (klon>1) ok_flux_surf=.false.
494 
495  ! Initilize debug IO
496  IF (debugindex .AND. mpi_size==1) THEN
497  ! initialize IOIPSL output
498  idayref = day_ini
499  CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
500  CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
501  DO i = 1, iim
502  zx_lon(i,1) = rlon(i+1)
503  zx_lon(i,jjm+1) = rlon(i+1)
504  ENDDO
505  CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
506  CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
507  1,iim,1,jjm+1, &
508  itau_phy,zjulian,dtime,nhoridbg,nidbg)
509  ! no vertical axis
510  cl_surf(1)='ter'
511  cl_surf(2)='lic'
512  cl_surf(3)='oce'
513  cl_surf(4)='sic'
514  DO nsrf=1,nbsrf
515  CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
516  jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
517  END DO
518 
519  CALL histend(nidbg)
520  CALL histsync(nidbg)
521 
522  END IF
523 
524  ENDIF
525 
526 !****************************************************************************************
527 ! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
528 ! instead of ORCHIDEE)
529  IF (qsol0>0.) THEN
530  print*,'WARNING : On impose qsol=',qsol0
531  qsol(:)=qsol0
532  ENDIF
533 !****************************************************************************************
534 
535 !****************************************************************************************
536 ! 2) Initialization to zero
537 ! Done for all local variables that will be compressed later
538 ! and argument with INTENT(OUT)
539 !****************************************************************************************
540  cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
541  ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
542  zv1 = 0.0 ; yqsurf = 0.0 ; yalb1 = 0.0 ; yalb2 = 0.0
543  yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
544  ysollw = 0.0 ; yrugos = 0.0 ; yu1 = 0.0
545  yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
546  ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
547  yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
548  yrugoro = 0.0 ; ywindsp = 0.0
549  d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
550  flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
551  d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.0 ; yqsol = 0.0
552  ytherm = 0.0 ; ytke=0.
553 
554  tke(:,:,is_ave)=0.
555  IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
556  zcoefh(:,:,:) = 0.0
557  zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
558  zcoefm(:,:,:) = 0.0
559  zcoefm(:,1,:) = 999999. !
560  ELSE
561  zcoefm(:,:,is_ave)=0.
562  zcoefh(:,:,is_ave)=0.
563  ENDIF
564  ytsoil = 999999.
565 
566  rh2m(:) = 0.
567  qsat2m(:) = 0.
568 !****************************************************************************************
569 ! 3) - Calculate pressure thickness of each layer
570 ! - Calculate the wind at first layer
571 ! - Mean calculations of albedo
572 ! - Calculate net radiance at sub-surface
573 !****************************************************************************************
574  DO k = 1, klev
575  DO i = 1, klon
576  delp(i,k) = paprs(i,k)-paprs(i,k+1)
577  ENDDO
578  ENDDO
579 
580 !****************************************************************************************
581 ! Test for rugos........ from physiq.. A la fin plutot???
582 !
583 !****************************************************************************************
584 
585  zxrugs(:) = 0.0
586  DO nsrf = 1, nbsrf
587  DO i = 1, klon
588  rugos(i,nsrf) = max(rugos(i,nsrf),0.000015)
589  zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
590  ENDDO
591  ENDDO
592 
593 ! Mean calculations of albedo
594 !
595 ! Albedo at sub-surface
596 ! * alb1 : albedo in visible SW interval
597 ! * alb2 : albedo in near infrared SW interval
598 ! * alb : mean albedo for whole SW interval
599 !
600 ! Mean albedo for grid point
601 ! * alb1_m : albedo in visible SW interval
602 ! * alb2_m : albedo in near infrared SW interval
603 ! * alb_m : mean albedo at whole SW interval
604 
605  alb1_m(:) = 0.0
606  alb2_m(:) = 0.0
607  DO nsrf = 1, nbsrf
608  DO i = 1, klon
609  alb1_m(i) = alb1_m(i) + alb1(i,nsrf) * pctsrf(i,nsrf)
610  alb2_m(i) = alb2_m(i) + alb2(i,nsrf) * pctsrf(i,nsrf)
611  ENDDO
612  ENDDO
613 
614 ! We here suppose the fraction f1 of incoming radiance of visible radiance
615 ! as a fraction of all shortwave radiance
616  f1 = 0.5
617 ! f1 = 1 ! put f1=1 to recreate old calculations
618 
619  DO nsrf = 1, nbsrf
620  DO i = 1, klon
621  alb(i,nsrf) = f1*alb1(i,nsrf) + (1-f1)*alb2(i,nsrf)
622  ENDDO
623  ENDDO
624 
625  DO i = 1, klon
626  alb_m(i) = f1*alb1_m(i) + (1-f1)*alb2_m(i)
627  END DO
628 
629 ! Calculation of mean temperature at surface grid points
630  ztsol(:) = 0.0
631  DO nsrf = 1, nbsrf
632  DO i = 1, klon
633  ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
634  ENDDO
635  ENDDO
636 
637 ! Linear distrubution on sub-surface of long- and shortwave net radiance
638  DO nsrf = 1, nbsrf
639  DO i = 1, klon
640  sollw(i,nsrf) = sollw_m(i) + 4.0*rsigma*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
641  solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
642  ENDDO
643  ENDDO
644 
645 
646 ! Downwelling longwave radiation at mean surface
647  lwdown_m(:) = 0.0
648  DO i = 1, klon
649  lwdown_m(i) = sollw_m(i) + rsigma*ztsol(i)**4
650  ENDDO
651 
652 !****************************************************************************************
653 ! 4) Loop over different surfaces
654 !
655 ! Only points containing a fraction of the sub surface will be threated.
656 !
657 !****************************************************************************************
658 
659  loop_nbsrf: DO nsrf = 1, nbsrf
660 
661 ! Search for index(ni) and size(knon) of domaine to treat
662  ni(:) = 0
663  knon = 0
664  DO i = 1, klon
665  IF (pctsrf(i,nsrf) > 0.) THEN
666  knon = knon + 1
667  ni(knon) = i
668  ENDIF
669  ENDDO
670 
671  ! write index, with IOIPSL
672  IF (debugindex .AND. mpi_size==1) THEN
673  tabindx(:)=0.
674  DO i=1,knon
675  tabindx(i)=REAL(i)
676  END DO
677  debugtab(:,:) = 0.
678  ndexbg(:) = 0
679  CALL gath2cpl(tabindx,debugtab,knon,ni)
680  CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
681  ENDIF
682 
683 !****************************************************************************************
684 ! 5) Compress variables
685 !
686 !****************************************************************************************
687 
688  DO j = 1, knon
689  i = ni(j)
690  ypct(j) = pctsrf(i,nsrf)
691  yts(j) = ts(i,nsrf)
692  ysnow(j) = snow(i,nsrf)
693  yqsurf(j) = qsurf(i,nsrf)
694  yalb(j) = alb(i,nsrf)
695  yalb1(j) = alb1(i,nsrf)
696  yalb2(j) = alb2(i,nsrf)
697  yrain_f(j) = rain_f(i)
698  ysnow_f(j) = snow_f(i)
699  yagesno(j) = agesno(i,nsrf)
700  yfder(j) = fder(i)
701  ysolsw(j) = solsw(i,nsrf)
702  ysollw(j) = sollw(i,nsrf)
703  yrugos(j) = rugos(i,nsrf)
704  yrugoro(j) = rugoro(i)
705  yu1(j) = u(i,1)
706  yv1(j) = v(i,1)
707  ypaprs(j,klev+1) = paprs(i,klev+1)
708  ywindsp(j) = sqrt(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )
709  END DO
710 
711  DO k = 1, klev
712  DO j = 1, knon
713  i = ni(j)
714  ypaprs(j,k) = paprs(i,k)
715  ypplay(j,k) = pplay(i,k)
716  ydelp(j,k) = delp(i,k)
717  ytke(j,k) = tke(i,k,nsrf)
718  yu(j,k) = u(i,k)
719  yv(j,k) = v(i,k)
720  yt(j,k) = t(i,k)
721  yq(j,k) = q(i,k)
722  ENDDO
723  ENDDO
724 
725  DO k = 1, nsoilmx
726  DO j = 1, knon
727  i = ni(j)
728  ytsoil(j,k) = ftsoil(i,k,nsrf)
729  END DO
730  END DO
731 
732  ! qsol(water height in soil) only for bucket continental model
733  IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
734  DO j = 1, knon
735  i = ni(j)
736  yqsol(j) = qsol(i)
737  END DO
738  ENDIF
739 
740 !****************************************************************************************
741 ! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
742 !
743 !****************************************************************************************
744 
745  CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
746  yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
747  yts, yqsurf, yrugos, &
748  ycdragm, ycdragh )
749 
750 !****************************************************************************************
751 ! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm.
752 !
753 !****************************************************************************************
754 
755  CALL coef_diff_turb(dtime, nsrf, knon, ni, &
756  ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
757  ycoefm, ycoefh, ytke)
758 
759  IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
760 ! In this case, coef_diff_turb is called for the Cd only
761  DO k = 2, klev
762  DO j = 1, knon
763  i = ni(j)
764  ycoefh(j,k) = zcoefh(i,k,nsrf)
765  ycoefm(j,k) = zcoefm(i,k,nsrf)
766  ENDDO
767  ENDDO
768  ENDIF
769 
770 !****************************************************************************************
771 !
772 ! 8) "La descente" - "The downhill"
773 !
774 ! climb_hq_down and climb_wind_down calculate the coefficients
775 ! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
776 ! Only the coefficients at surface for H and Q are returned.
777 !
778 !****************************************************************************************
779 
780 ! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
781  CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
782  ydelp, yt, yq, dtime, &
783  acoefh, acoefq, bcoefh, bcoefq)
784 
785 ! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
786  CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
787  acoefu, acoefv, bcoefu, bcoefv)
788 
789 
790 !****************************************************************************************
791 ! 9) Small calculations
792 !
793 !****************************************************************************************
794 
795 ! - Reference pressure is given the values at surface level
796  ypsref(:) = ypaprs(:,1)
797 
798 ! - CO2 field on 2D grid to be sent to ORCHIDEE
799 ! Transform to compressed field
800  IF (carbon_cycle_cpl) THEN
801  DO i=1,knon
802  r_co2_ppm(i) = co2_send(ni(i))
803  END DO
804  ELSE
805  r_co2_ppm(:) = co2_ppm ! Constant field
806  END IF
807 
808 !****************************************************************************************
809 !
810 ! Calulate t2m and q2m for the case of calculation at land grid points
811 ! t2m and q2m are needed as input to ORCHIDEE
812 !
813 !****************************************************************************************
814  IF (nsrf == is_ter) THEN
815 
816  DO i = 1, knon
817  zgeo1(i) = rd * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
818  * (ypaprs(i,1)-ypplay(i,1))
819  END DO
820 
821  ! Calculate the temperature et relative humidity at 2m and the wind at 10m
822  CALL stdlevvar(klon, knon, is_ter, zxli, &
823  yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
824  yts, yqsurf, yrugos, ypaprs(:,1), ypplay(:,1), &
825  yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
826 
827  END IF
828 
829 !****************************************************************************************
830 !
831 ! 10) Switch selon current surface
832 ! It is necessary to start with the continental surfaces because the ocean
833 ! needs their run-off.
834 !
835 !****************************************************************************************
836  SELECT CASE(nsrf)
837 
838  CASE(is_ter)
839  ! ylwdown : to be removed, calculation is now done at land surface in surf_land
840  ylwdown(:)=0.0
841  DO i=1,knon
842  ylwdown(i)=lwdown_m(ni(i))
843  END DO
844  CALL surf_land(itap, dtime, date0, jour, knon, ni,&
845  rlon, rlat, &
846  debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
847  yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
848  acoefh, acoefq, bcoefh, bcoefq, &
849  acoefu, acoefv, bcoefu, bcoefv, &
850  ypsref, yu1, yv1, yrugoro, pctsrf, &
851  ylwdown, yq2m, yt2m, &
852  ysnow, yqsol, yagesno, ytsoil, &
853  yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
854  yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
855  y_flux_u1, y_flux_v1 )
856 
857 
858  CASE(is_lic)
859  CALL surf_landice(itap, dtime, knon, ni, &
860  ysolsw, ysollw, yts, ypplay(:,1), &
861  ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
862  acoefh, acoefq, bcoefh, bcoefq, &
863  acoefu, acoefv, bcoefu, bcoefv, &
864  ypsref, yu1, yv1, yrugoro, pctsrf, &
865  ysnow, yqsurf, yqsol, yagesno, &
866  ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
867  ytsurf_new, y_dflux_t, y_dflux_q, &
868  y_flux_u1, y_flux_v1)
869 
870  CASE(is_oce)
871  CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
872  yrugos, ywindsp, rmu0, yfder, yts, &
873  itap, dtime, jour, knon, ni, &
874  ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
875  acoefh, acoefq, bcoefh, bcoefq, &
876  acoefu, acoefv, bcoefu, bcoefv, &
877  ypsref, yu1, yv1, yrugoro, pctsrf, &
878  ysnow, yqsurf, yagesno, &
879  yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
880  ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
881  y_flux_u1, y_flux_v1)
882 
883  CASE(is_sic)
884  CALL surf_seaice( &
885  rlon, rlat, ysolsw, ysollw, yalb1, yfder, &
886  itap, dtime, jour, knon, ni, &
887  lafin, &
888  yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
889  acoefh, acoefq, bcoefh, bcoefq, &
890  acoefu, acoefv, bcoefu, bcoefv, &
891  ypsref, yu1, yv1, yrugoro, pctsrf, &
892  ysnow, yqsurf, yqsol, yagesno, ytsoil, &
893  yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
894  ytsurf_new, y_dflux_t, y_dflux_q, &
895  y_flux_u1, y_flux_v1)
896 
897 
898  CASE default
899  WRITE(lunout,*) 'Surface index = ', nsrf
900  abort_message = 'Surface index not valid'
901  CALL abort_gcm(modname,abort_message,1)
902  END SELECT
903 
904 
905 !****************************************************************************************
906 ! 11) - Calcul the increment of surface temperature
907 !
908 !****************************************************************************************
909  y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
910 
911 !****************************************************************************************
912 !
913 ! 12) "La remontee" - "The uphill"
914 !
915 ! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
916 ! for X=H, Q, U and V, for all vertical levels.
917 !
918 !****************************************************************************************
919 ! H and Q
920  IF (ok_flux_surf) THEN
921  print *,'pbl_surface: fsens flat RLVTT=',fsens,flat,rlvtt
922  y_flux_t1(:) = fsens
923  y_flux_q1(:) = flat/rlvtt
924  yfluxlat(:) = flat
925 
926  kech_h(:) = ycdragh(:) * (1.0+sqrt(yu(:,1)**2+yv(:,1)**2)) * &
927  ypplay(:,1)/(rd*yt(:,1))
928  ytoto(:)=(1./rcpd)*(acoefh(:)+bcoefh(:)*y_flux_t1(:)*dtime)
929  ytsurf_new(:)=ytoto(:)-y_flux_t1(:)/(kech_h(:)*rcpd)
930  y_d_ts(:) = ytsurf_new(:) - yts(:)
931 
932  ELSE
933  y_flux_t1(:) = yfluxsens(:)
934  y_flux_q1(:) = -yevap(:)
935  ENDIF
936 
937  CALL climb_hq_up(knon, dtime, yt, yq, &
938  y_flux_q1, y_flux_t1, ypaprs, ypplay, &
939  y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
940 
941 
942  CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
943  y_flux_u, y_flux_v, y_d_u, y_d_v)
944 
945 
946  y_d_t_diss(:,:)=0.
947  IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
948  CALL yamada_c(knon,dtime,ypaprs,ypplay &
949  & ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
950  & ,iflag_pbl,nsrf)
951  ENDIF
952 ! print*,'yamada_c OK'
953 
954  DO j = 1, knon
955  y_dflux_t(j) = y_dflux_t(j) * ypct(j)
956  y_dflux_q(j) = y_dflux_q(j) * ypct(j)
957  ENDDO
958 
959 !****************************************************************************************
960 ! 13) Transform variables for output format :
961 ! - Decompress
962 ! - Multiply with pourcentage of current surface
963 ! - Cumulate in global variable
964 !
965 !****************************************************************************************
966 
967  DO k = 1, klev
968  DO j = 1, knon
969  i = ni(j)
970  y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)
971  y_d_t(j,k) = y_d_t(j,k) * ypct(j)
972  y_d_q(j,k) = y_d_q(j,k) * ypct(j)
973  y_d_u(j,k) = y_d_u(j,k) * ypct(j)
974  y_d_v(j,k) = y_d_v(j,k) * ypct(j)
975 
976  flux_t(i,k,nsrf) = y_flux_t(j,k)
977  flux_q(i,k,nsrf) = y_flux_q(j,k)
978  flux_u(i,k,nsrf) = y_flux_u(j,k)
979  flux_v(i,k,nsrf) = y_flux_v(j,k)
980 
981 
982  ENDDO
983  ENDDO
984 
985 ! print*,'Dans pbl OK1'
986 
987  evap(:,nsrf) = - flux_q(:,1,nsrf)
988 
989  alb1(:, nsrf) = 0.
990  alb2(:, nsrf) = 0.
991  snow(:, nsrf) = 0.
992  qsurf(:, nsrf) = 0.
993  rugos(:, nsrf) = 0.
994  fluxlat(:,nsrf) = 0.
995  DO j = 1, knon
996  i = ni(j)
997  d_ts(i,nsrf) = y_d_ts(j)
998  alb1(i,nsrf) = yalb1_new(j)
999  alb2(i,nsrf) = yalb2_new(j)
1000  snow(i,nsrf) = ysnow(j)
1001  qsurf(i,nsrf) = yqsurf(j)
1002  rugos(i,nsrf) = yz0_new(j)
1003  fluxlat(i,nsrf) = yfluxlat(j)
1004  agesno(i,nsrf) = yagesno(j)
1005  cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
1006  cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
1007  dflux_t(i) = dflux_t(i) + y_dflux_t(j)
1008  dflux_q(i) = dflux_q(i) + y_dflux_q(j)
1009  END DO
1010 
1011 ! print*,'Dans pbl OK2'
1012 
1013  DO k = 2, klev
1014  DO j = 1, knon
1015  i = ni(j)
1016  tke(i,k,nsrf) = ytke(j,k)
1017  zcoefh(i,k,nsrf) = ycoefh(j,k)
1018  zcoefm(i,k,nsrf) = ycoefm(j,k)
1019  tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
1020  zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
1021  zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
1022  END DO
1023  END DO
1024 
1025 ! print*,'Dans pbl OK3'
1026 
1027  IF ( nsrf .EQ. is_ter ) THEN
1028  DO j = 1, knon
1029  i = ni(j)
1030  qsol(i) = yqsol(j)
1031  END DO
1032  END IF
1033 
1034  ftsoil(:,:,nsrf) = 0.
1035  DO k = 1, nsoilmx
1036  DO j = 1, knon
1037  i = ni(j)
1038  ftsoil(i, k, nsrf) = ytsoil(j,k)
1039  END DO
1040  END DO
1041 
1042 
1043  DO k = 1, klev
1044  DO j = 1, knon
1045  i = ni(j)
1046  d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
1047  d_t(i,k) = d_t(i,k) + y_d_t(j,k)
1048  d_q(i,k) = d_q(i,k) + y_d_q(j,k)
1049  d_u(i,k) = d_u(i,k) + y_d_u(j,k)
1050  d_v(i,k) = d_v(i,k) + y_d_v(j,k)
1051  END DO
1052  END DO
1053 
1054 ! print*,'Dans pbl OK4'
1055 
1056 !****************************************************************************************
1057 ! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
1058 ! Call HBTM
1059 !
1060 !****************************************************************************************
1061  t2m(:,nsrf) = 0.
1062  q2m(:,nsrf) = 0.
1063  ustar(:,nsrf) = 0.
1064  u10m(:,nsrf) = 0.
1065  v10m(:,nsrf) = 0.
1066  pblh(:,nsrf) = 0. ! Hauteur de couche limite
1067  plcl(:,nsrf) = 0. ! Niveau de condensation de la CLA
1068  capcl(:,nsrf) = 0. ! CAPE de couche limite
1069  oliqcl(:,nsrf) = 0. ! eau_liqu integree de couche limite
1070  cteicl(:,nsrf) = 0. ! cloud top instab. crit. couche limite
1071  pblt(:,nsrf) = 0. ! T a la Hauteur de couche limite
1072  therm(:,nsrf) = 0.
1073  trmb1(:,nsrf) = 0. ! deep_cape
1074  trmb2(:,nsrf) = 0. ! inhibition
1075  trmb3(:,nsrf) = 0. ! Point Omega
1076 
1077 #undef T2m
1078 #define T2m
1079 #ifdef T2m
1080 ! Calculations of diagnostic t,q at 2m and u, v at 10m
1081 
1082 ! print*,'Dans pbl OK41'
1083 ! print*,'tair1,yt(:,1),y_d_t(:,1)'
1084 ! print*, tair1,yt(:,1),y_d_t(:,1)
1085  DO j=1, knon
1086  i = ni(j)
1087  uzon(j) = yu(j,1) + y_d_u(j,1)
1088  vmer(j) = yv(j,1) + y_d_v(j,1)
1089  tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
1090  qair1(j) = yq(j,1) + y_d_q(j,1)
1091  zgeo1(j) = rd * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
1092  * (ypaprs(j,1)-ypplay(j,1))
1093  tairsol(j) = yts(j) + y_d_ts(j)
1094  rugo1(j) = yrugos(j)
1095  IF(nsrf.EQ.is_oce) THEN
1096  rugo1(j) = rugos(i,nsrf)
1097  ENDIF
1098  psfce(j)=ypaprs(j,1)
1099  patm(j)=ypplay(j,1)
1100  qairsol(j) = yqsurf(j)
1101  END DO
1102 
1103 ! print*,'Dans pbl OK42A'
1104 ! print*,'tair1,yt(:,1),y_d_t(:,1)'
1105 ! print*, tair1,yt(:,1),y_d_t(:,1)
1106 
1107 ! Calculate the temperature et relative humidity at 2m and the wind at 10m
1108  CALL stdlevvar(klon, knon, nsrf, zxli, &
1109  uzon, vmer, tair1, qair1, zgeo1, &
1110  tairsol, qairsol, rugo1, psfce, patm, &
1111  yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1112 ! print*,'Dans pbl OK42B'
1113 
1114  DO j=1, knon
1115  i = ni(j)
1116  t2m(i,nsrf)=yt2m(j)
1117  q2m(i,nsrf)=yq2m(j)
1118 
1119  ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
1120  ustar(i,nsrf)=yustar(j)
1121  u10m(i,nsrf)=(yu10m(j) * uzon(j))/sqrt(uzon(j)**2+vmer(j)**2)
1122  v10m(i,nsrf)=(yu10m(j) * vmer(j))/sqrt(uzon(j)**2+vmer(j)**2)
1123 
1124  END DO
1125 
1126 ! print*,'Dans pbl OK43'
1127 !IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
1128 !IM Ajoute dependance type surface
1129  IF (thermcep) THEN
1130  DO j = 1, knon
1131  i=ni(j)
1132  zdelta1 = max(0.,sign(1., rtt-yt2m(j) ))
1133  zx_qs1 = r2es * foeew(yt2m(j),zdelta1)/paprs(i,1)
1134  zx_qs1 = min(0.5,zx_qs1)
1135  zcor1 = 1./(1.-retv*zx_qs1)
1136  zx_qs1 = zx_qs1*zcor1
1137 
1138  rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
1139  qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
1140  END DO
1141  END IF
1142 
1143 ! print*,'OK pbl 5'
1144  CALL hbtm(knon, ypaprs, ypplay, &
1145  yt2m,yt10m,yq2m,yq10m,yustar, &
1146  y_flux_t,y_flux_q,yu,yv,yt,yq, &
1147  ypblh,ycapcl,yoliqcl,ycteicl,ypblt, &
1148  ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
1149 
1150  DO j=1, knon
1151  i = ni(j)
1152  pblh(i,nsrf) = ypblh(j)
1153  plcl(i,nsrf) = ylcl(j)
1154  capcl(i,nsrf) = ycapcl(j)
1155  oliqcl(i,nsrf) = yoliqcl(j)
1156  cteicl(i,nsrf) = ycteicl(j)
1157  pblt(i,nsrf) = ypblt(j)
1158  therm(i,nsrf) = ytherm(j)
1159  trmb1(i,nsrf) = ytrmb1(j)
1160  trmb2(i,nsrf) = ytrmb2(j)
1161  trmb3(i,nsrf) = ytrmb3(j)
1162  END DO
1163 
1164 ! print*,'OK pbl 6'
1165 #else
1166 ! T2m not defined
1167 ! No calculation
1168  print*,' Warning !!! No T2m calculation. Output is set to zero.'
1169 #endif
1170 
1171 !****************************************************************************************
1172 ! 15) End of loop over different surfaces
1173 !
1174 !****************************************************************************************
1175  END DO loop_nbsrf
1176 
1177 !****************************************************************************************
1178 ! 16) Calculate the mean value over all sub-surfaces for som variables
1179 !
1180 !****************************************************************************************
1181 
1182 ! print*,'OK pbl 7'
1183  zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1184  zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1185  DO nsrf = 1, nbsrf
1186  DO k = 1, klev
1187  DO i = 1, klon
1188  zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
1189  zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
1190  zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
1191  zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
1192  END DO
1193  END DO
1194  END DO
1195 
1196 ! print*,'OK pbl 8'
1197  DO i = 1, klon
1198  zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1199  zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1200  fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1201  ENDDO
1202 
1203 !
1204 ! Incrementer la temperature du sol
1205 !
1206  zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
1207  zt2m(:) = 0.0 ; zq2m(:) = 0.0
1208  zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
1209  s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
1210  s_capcl(:) = 0.0 ; s_oliqcl(:) = 0.0
1211  s_cteicl(:) = 0.0; s_pblt(:) = 0.0
1212  s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1213  s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1214 
1215 ! print*,'OK pbl 9'
1216 
1217  DO nsrf = 1, nbsrf
1218  DO i = 1, klon
1219  ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1220 
1221  wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1222  + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1223  wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1224  pctsrf(i,nsrf)
1225 
1226  zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
1227  zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
1228 
1229  zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
1230  zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
1231  zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
1232  zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
1233  zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
1234 
1235  s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
1236  s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
1237  s_capcl(i) = s_capcl(i) + capcl(i,nsrf) * pctsrf(i,nsrf)
1238  s_oliqcl(i) = s_oliqcl(i) + oliqcl(i,nsrf)* pctsrf(i,nsrf)
1239  s_cteicl(i) = s_cteicl(i) + cteicl(i,nsrf)* pctsrf(i,nsrf)
1240  s_pblt(i) = s_pblt(i) + pblt(i,nsrf) * pctsrf(i,nsrf)
1241  s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
1242  s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
1243  s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
1244  s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
1245  END DO
1246  END DO
1247 ! print*,'OK pbl 10'
1248 
1249  IF (check) THEN
1250  amn=min(ts(1,is_ter),1000.)
1251  amx=max(ts(1,is_ter),-1000.)
1252  DO i=2, klon
1253  amn=min(ts(i,is_ter),amn)
1254  amx=max(ts(i,is_ter),amx)
1255  ENDDO
1256  print*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1257  ENDIF
1258 
1259 !jg ?
1260 !!$!
1261 !!$! If a sub-surface does not exsist for a grid point, the mean value for all
1262 !!$! sub-surfaces is distributed.
1263 !!$!
1264 !!$ DO nsrf = 1, nbsrf
1265 !!$ DO i = 1, klon
1266 !!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1267 !!$ ts(i,nsrf) = zxtsol(i)
1268 !!$ t2m(i,nsrf) = zt2m(i)
1269 !!$ q2m(i,nsrf) = zq2m(i)
1270 !!$ u10m(i,nsrf) = zu10m(i)
1271 !!$ v10m(i,nsrf) = zv10m(i)
1272 !!$
1273 !!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1274 !!$ pblh(i,nsrf) = s_pblh(i)
1275 !!$ plcl(i,nsrf) = s_plcl(i)
1276 !!$ capCL(i,nsrf) = s_capCL(i)
1277 !!$ oliqCL(i,nsrf) = s_oliqCL(i)
1278 !!$ cteiCL(i,nsrf) = s_cteiCL(i)
1279 !!$ pblT(i,nsrf) = s_pblT(i)
1280 !!$ therm(i,nsrf) = s_therm(i)
1281 !!$ trmb1(i,nsrf) = s_trmb1(i)
1282 !!$ trmb2(i,nsrf) = s_trmb2(i)
1283 !!$ trmb3(i,nsrf) = s_trmb3(i)
1284 !!$ ENDIF
1285 !!$ ENDDO
1286 !!$ ENDDO
1287 
1288 
1289  DO i = 1, klon
1290  fder(i) = - 4.0*rsigma*zxtsol(i)**3
1291  ENDDO
1292 
1293  zxqsurf(:) = 0.0
1294  zxsnow(:) = 0.0
1295  DO nsrf = 1, nbsrf
1296  DO i = 1, klon
1297  zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
1298  zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
1299  END DO
1300  END DO
1301 
1302 ! Premier niveau de vent sortie dans physiq.F
1303  zu1(:) = u(:,1)
1304  zv1(:) = v(:,1)
1305 
1306 ! Some of the module declared variables are returned for printing in physiq.F
1307  qsol_d(:) = qsol(:)
1308  evap_d(:,:) = evap(:,:)
1309  rugos_d(:,:) = rugos(:,:)
1310  agesno_d(:,:) = agesno(:,:)
1311 
1312 
1313  END SUBROUTINE pbl_surface
1314 !
1315 !****************************************************************************************
1316 !
1317  SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1318  evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1319 
1320  include "indicesol.h"
1321  include "dimsoil.h"
1322 
1323 ! Ouput variables
1324 !****************************************************************************************
1325  REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst
1326  REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
1327  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
1328  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
1329  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst
1330  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst
1331  REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst
1332  REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1333 
1334 
1335 !****************************************************************************************
1336 ! Return module variables for writing to restart file
1337 !
1338 !****************************************************************************************
1339  qsol_rst(:) = qsol(:)
1340  fder_rst(:) = fder(:)
1341  snow_rst(:,:) = snow(:,:)
1342  qsurf_rst(:,:) = qsurf(:,:)
1343  evap_rst(:,:) = evap(:,:)
1344  rugos_rst(:,:) = rugos(:,:)
1345  agesno_rst(:,:) = agesno(:,:)
1346  ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1347 
1348 !****************************************************************************************
1349 ! Deallocate module variables
1350 !
1351 !****************************************************************************************
1352 ! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1353  IF (ALLOCATED(qsol)) DEALLOCATE(qsol)
1354  IF (ALLOCATED(fder)) DEALLOCATE(fder)
1355  IF (ALLOCATED(snow)) DEALLOCATE(snow)
1356  IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
1357  IF (ALLOCATED(evap)) DEALLOCATE(evap)
1358  IF (ALLOCATED(rugos)) DEALLOCATE(rugos)
1359  IF (ALLOCATED(agesno)) DEALLOCATE(agesno)
1360  IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
1361 
1362  END SUBROUTINE pbl_surface_final
1363 !
1364 !****************************************************************************************
1365 !
1366  SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, ustar, u10m, v10m, tke)
1367 
1368  ! Give default values where new fraction has appread
1369 
1370  include "indicesol.h"
1371  include "dimsoil.h"
1372  include "clesphys.h"
1373  include "compbl.h"
1374 
1375 ! Input variables
1376 !****************************************************************************************
1377  INTEGER, INTENT(IN) :: itime
1378  REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
1379 
1380 ! InOutput variables
1381 !****************************************************************************************
1382  REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
1383  REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb1, alb2
1384  REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
1385  REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
1386 
1387 ! Local variables
1388 !****************************************************************************************
1389  INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
1390  CHARACTER(len=80) :: abort_message
1391  CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
1392  INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
1393 !
1394 ! All at once !!
1395 !****************************************************************************************
1396 
1397  DO nsrf = 1, nbsrf
1398  ! First decide complement sub-surfaces
1399  SELECT CASE (nsrf)
1400  CASE(is_oce)
1401  nsrf_comp1=is_sic
1402  nsrf_comp2=is_ter
1403  nsrf_comp3=is_lic
1404  CASE(is_sic)
1405  nsrf_comp1=is_oce
1406  nsrf_comp2=is_ter
1407  nsrf_comp3=is_lic
1408  CASE(is_ter)
1409  nsrf_comp1=is_lic
1410  nsrf_comp2=is_oce
1411  nsrf_comp3=is_sic
1412  CASE(is_lic)
1413  nsrf_comp1=is_ter
1414  nsrf_comp2=is_oce
1415  nsrf_comp3=is_sic
1416  END SELECT
1417 
1418  ! Initialize all new fractions
1419  DO i=1, klon
1420  IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
1421 
1422  IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
1423  ! Use the complement sub-surface, keeping the continents unchanged
1424  qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
1425  evap(i,nsrf) = evap(i,nsrf_comp1)
1426  rugos(i,nsrf) = rugos(i,nsrf_comp1)
1427  tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
1428  alb1(i,nsrf) = alb1(i,nsrf_comp1)
1429  alb2(i,nsrf) = alb2(i,nsrf_comp1)
1430  ustar(i,nsrf) = ustar(i,nsrf_comp1)
1431  u10m(i,nsrf) = u10m(i,nsrf_comp1)
1432  v10m(i,nsrf) = v10m(i,nsrf_comp1)
1433  if (iflag_pbl > 1) then
1434  tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
1435  endif
1436  mfois(nsrf) = mfois(nsrf) + 1
1437  ELSE
1438  ! The continents have changed. The new fraction receives the mean sum of the existent fractions
1439  qsurf(i,nsrf) = qsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
1440  evap(i,nsrf) = evap(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1441  rugos(i,nsrf) = rugos(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
1442  tsurf(i,nsrf) = tsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
1443  alb1(i,nsrf) = alb1(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1444  alb2(i,nsrf) = alb2(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1445  ustar(i,nsrf) = ustar(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1446  u10m(i,nsrf) = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1447  v10m(i,nsrf) = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1448  if (iflag_pbl > 1) then
1449  tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
1450  endif
1451 
1452  ! Security abort. This option has never been tested. To test, comment the following line.
1453 ! abort_message='The fraction of the continents have changed!'
1454 ! CALL abort_gcm(modname,abort_message,1)
1455  nfois(nsrf) = nfois(nsrf) + 1
1456  END IF
1457  snow(i,nsrf) = 0.
1458  agesno(i,nsrf) = 0.
1459  ftsoil(i,:,nsrf) = tsurf(i,nsrf)
1460  ELSE
1461  pfois(nsrf) = pfois(nsrf)+ 1
1462  END IF
1463  END DO
1464 
1465  END DO
1466 
1467  END SUBROUTINE pbl_surface_newfrac
1468 
1469 !
1470 !****************************************************************************************
1471 !
1472 
1473 END MODULE pbl_surface_mod