source: LMDZ5/branches/LMDZ5_SPLA/libf/phylmd/pbl_surface_mod.F90 @ 2404

Last change on this file since 2404 was 1932, checked in by lguez, 11 years ago

Declaration of variables appearing in array bounds must be specified
before the declaration of the array.

Variables used in FCTTRE.h are declared in YOMCST.h and YOETHF.h so
YOMCST.h and YOETHF.h must be included before FCTTRE.h.

(There were compilation errors with gfortran and debugging options.)

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