source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/pbl_surface_mod.F90 @ 3655

Last change on this file since 3655 was 1278, checked in by Laurent Fairhead, 15 years ago

Pour le 1D MPL

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