source: LMDZ5/branches/LF-private/libf/phylmd/pbl_surface_mod.F90 @ 2261

Last change on this file since 2261 was 1887, checked in by Laurent Fairhead, 11 years ago

Common ought to be placed in an include file

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