source: LMDZ4/trunk/libf/phylmd/pbl_surface_mod.F90 @ 1357

Last change on this file since 1357 was 1282, checked in by Ehouarn Millour, 15 years ago

Bug fix: variables cannot be both 'save' and 'common'

(this error is tolerated by most compilers;

but not xlf)

EM

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 58.2 KB
Line 
1!
2! $Id: pbl_surface_mod.F90 1282 2009-12-16 11:23:05Z emillour $
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 ! initialized during first_call below
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       ! Initialize ok_flux_surf (for 1D model)
485       ok_flux_surf=.FALSE.
486       
487       ! Initilize debug IO
488       IF (debugindex .AND. mpi_size==1) THEN
489          ! initialize IOIPSL output
490          idayref = day_ini
491          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
492          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
493          DO i = 1, iim
494             zx_lon(i,1) = rlon(i+1)
495             zx_lon(i,jjm+1) = rlon(i+1)
496          ENDDO
497          CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
498          CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
499               1,iim,1,jjm+1, &
500               itau_phy,zjulian,dtime,nhoridbg,nidbg)
501          ! no vertical axis
502          cl_surf(1)='ter'
503          cl_surf(2)='lic'
504          cl_surf(3)='oce'
505          cl_surf(4)='sic'
506          DO nsrf=1,nbsrf
507             CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
508                  jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
509          END DO
510
511          CALL histend(nidbg)
512          CALL histsync(nidbg)
513
514       END IF
515       
516    ENDIF
517         
518!****************************************************************************************
519! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
520! instead of ORCHIDEE)
521    IF (qsol0>0.) THEN
522      PRINT*,'WARNING : On impose qsol=',qsol0
523      qsol(:)=qsol0
524    ENDIF
525!****************************************************************************************
526
527!****************************************************************************************
528! 2) Initialization to zero
529!    Done for all local variables that will be compressed later
530!    and argument with INTENT(OUT)
531!****************************************************************************************
532    cdragh = 0.0  ; cdragm = 0.0     ; dflux_t = 0.0   ; dflux_q = 0.0
533    ypct = 0.0    ; yts = 0.0        ; ysnow = 0.0
534    zv1 = 0.0     ; yqsurf = 0.0     ; yalb1 = 0.0     ; yalb2 = 0.0   
535    yrain_f = 0.0 ; ysnow_f = 0.0    ; yfder = 0.0     ; ysolsw = 0.0   
536    ysollw = 0.0  ; yrugos = 0.0     ; yu1 = 0.0   
537    yv1 = 0.0     ; ypaprs = 0.0     ; ypplay = 0.0
538    ydelp = 0.0   ; yu = 0.0         ; yv = 0.0        ; yt = 0.0         
539    yq = 0.0      ; y_dflux_t = 0.0  ; y_dflux_q = 0.0
540    yrugoro = 0.0 ; ywindsp = 0.0   
541    d_ts = 0.0    ; yfluxlat=0.0     ; flux_t = 0.0    ; flux_q = 0.0     
542    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
543    d_u = 0.0     ; d_v = 0.0        ; yqsol = 0.0   
544    ytherm = 0.0  ; ytke=0.
545   
546    zcoefh(:,:) = 0.0
547    zcoefh(:,1) = 999999. ! zcoefh(:,k=1) should never be used
548    ytsoil = 999999.
549
550    rh2m(:)        = 0.
551    qsat2m(:)      = 0.
552!****************************************************************************************
553! 3) - Calculate pressure thickness of each layer
554!    - Calculate the wind at first layer
555!    - Mean calculations of albedo
556!    - Calculate net radiance at sub-surface
557!****************************************************************************************
558    DO k = 1, klev
559       DO i = 1, klon
560          delp(i,k) = paprs(i,k)-paprs(i,k+1)
561       ENDDO
562    ENDDO
563
564!****************************************************************************************
565! Test for rugos........ from physiq.. A la fin plutot???
566!
567!****************************************************************************************
568
569    zxrugs(:) = 0.0
570    DO nsrf = 1, nbsrf
571       DO i = 1, klon
572          rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
573          zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
574       ENDDO
575    ENDDO
576
577! Mean calculations of albedo
578!
579! Albedo at sub-surface
580! * alb1 : albedo in visible SW interval
581! * alb2 : albedo in near infrared SW interval
582! * alb  : mean albedo for whole SW interval
583!
584! Mean albedo for grid point
585! * alb1_m : albedo in visible SW interval
586! * alb2_m : albedo in near infrared SW interval
587! * alb_m  : mean albedo at whole SW interval
588
589    alb1_m(:) = 0.0
590    alb2_m(:) = 0.0
591    DO nsrf = 1, nbsrf
592       DO i = 1, klon
593          alb1_m(i) = alb1_m(i) + alb1(i,nsrf) * pctsrf(i,nsrf)
594          alb2_m(i) = alb2_m(i) + alb2(i,nsrf) * pctsrf(i,nsrf)
595       ENDDO
596    ENDDO
597
598! We here suppose the fraction f1 of incoming radiance of visible radiance
599! as a fraction of all shortwave radiance
600    f1 = 0.5
601!    f1 = 1    ! put f1=1 to recreate old calculations
602
603    DO nsrf = 1, nbsrf
604       DO i = 1, klon
605          alb(i,nsrf) = f1*alb1(i,nsrf) + (1-f1)*alb2(i,nsrf)
606       ENDDO
607    ENDDO
608
609    DO i = 1, klon
610       alb_m(i) = f1*alb1_m(i) + (1-f1)*alb2_m(i)
611    END DO
612
613! Calculation of mean temperature at surface grid points
614    ztsol(:) = 0.0
615    DO nsrf = 1, nbsrf
616       DO i = 1, klon
617          ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
618       ENDDO
619    ENDDO
620
621! Linear distrubution on sub-surface of long- and shortwave net radiance
622    DO nsrf = 1, nbsrf
623       DO i = 1, klon
624          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
625          solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
626       ENDDO
627    ENDDO
628
629
630! Downwelling longwave radiation at mean surface
631    lwdown_m(:) = 0.0
632    DO i = 1, klon
633       lwdown_m(i) = sollw_m(i) + RSIGMA*ztsol(i)**4
634    ENDDO
635
636!****************************************************************************************
637! 4) Loop over different surfaces
638!
639! Only points containing a fraction of the sub surface will be threated.
640!
641!****************************************************************************************
642   
643    loop_nbsrf: DO nsrf = 1, nbsrf
644
645! Search for index(ni) and size(knon) of domaine to treat
646       ni(:) = 0
647       knon  = 0
648       DO i = 1, klon
649          IF (pctsrf(i,nsrf) > 0.) THEN
650             knon = knon + 1
651             ni(knon) = i
652          ENDIF
653       ENDDO
654
655       ! write index, with IOIPSL
656       IF (debugindex .AND. mpi_size==1) THEN
657          tabindx(:)=0.
658          DO i=1,knon
659             tabindx(i)=FLOAT(i)
660          END DO
661          debugtab(:,:) = 0.
662          ndexbg(:) = 0
663          CALL gath2cpl(tabindx,debugtab,knon,ni)
664          CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
665       ENDIF
666       
667!****************************************************************************************
668! 5) Compress variables
669!
670!****************************************************************************************
671
672       DO j = 1, knon
673          i = ni(j)
674          ypct(j)    = pctsrf(i,nsrf)
675          yts(j)     = ts(i,nsrf)
676          ysnow(j)   = snow(i,nsrf)
677          yqsurf(j)  = qsurf(i,nsrf)
678          yalb(j)    = alb(i,nsrf)
679          yalb1(j)   = alb1(i,nsrf)
680          yalb2(j)   = alb2(i,nsrf)
681          yrain_f(j) = rain_f(i)
682          ysnow_f(j) = snow_f(i)
683          yagesno(j) = agesno(i,nsrf)
684          yfder(j)   = fder(i)
685          ysolsw(j)  = solsw(i,nsrf)
686          ysollw(j)  = sollw(i,nsrf)
687          yrugos(j)  = rugos(i,nsrf)
688          yrugoro(j) = rugoro(i)
689          yu1(j)     = u(i,1)
690          yv1(j)     = v(i,1)
691          ypaprs(j,klev+1) = paprs(i,klev+1)
692          ywindsp(j) = SQRT(u10m(i,nsrf)**2 + v10m(i,nsrf)**2 )
693       END DO
694
695       DO k = 1, klev
696          DO j = 1, knon
697             i = ni(j)
698             ypaprs(j,k) = paprs(i,k)
699             ypplay(j,k) = pplay(i,k)
700             ydelp(j,k)  = delp(i,k)
701             ytke(j,k)   = tke(i,k,nsrf)
702             yu(j,k) = u(i,k)
703             yv(j,k) = v(i,k)
704             yt(j,k) = t(i,k)
705             yq(j,k) = q(i,k)
706          ENDDO
707       ENDDO
708       
709       DO k = 1, nsoilmx
710          DO j = 1, knon
711             i = ni(j)
712             ytsoil(j,k) = ftsoil(i,k,nsrf)
713          END DO
714       END DO
715       
716       ! qsol(water height in soil) only for bucket continental model
717       IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
718          DO j = 1, knon
719             i = ni(j)
720             yqsol(j) = qsol(i)
721          END DO
722       ENDIF
723       
724!****************************************************************************************
725! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
726!
727!****************************************************************************************
728
729       CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
730            yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
731            yts, yqsurf, yrugos, &
732            ycdragm, ycdragh )
733
734!****************************************************************************************
735! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm.
736!
737!****************************************************************************************
738
739       CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
740            ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
741            ycoefm, ycoefh, ytke)
742       
743!****************************************************************************************
744!
745! 8) "La descente" - "The downhill"
746
747!  climb_hq_down and climb_wind_down calculate the coefficients
748!  Ccoef_X et Dcoef_X for X=[H, Q, U, V].
749!  Only the coefficients at surface for H and Q are returned.
750!
751!****************************************************************************************
752
753! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
754       CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
755            ydelp, yt, yq, dtime, &
756            AcoefH, AcoefQ, BcoefH, BcoefQ)
757
758! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
759       CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
760            AcoefU, AcoefV, BcoefU, BcoefV)
761     
762
763!****************************************************************************************
764! 9) Small calculations
765!
766!****************************************************************************************
767
768! - Reference pressure is given the values at surface level         
769       ypsref(:) = ypaprs(:,1) 
770
771! - CO2 field on 2D grid to be sent to ORCHIDEE
772!   Transform to compressed field
773       IF (carbon_cycle_cpl) THEN
774          DO i=1,knon
775             r_co2_ppm(i) = co2_send(ni(i))
776          END DO
777       ELSE
778          r_co2_ppm(:) = co2_ppm     ! Constant field
779       END IF
780
781!****************************************************************************************
782!
783! Calulate t2m and q2m for the case of calculation at land grid points
784! t2m and q2m are needed as input to ORCHIDEE
785!
786!****************************************************************************************
787       IF (nsrf == is_ter) THEN
788
789          DO i = 1, knon
790             zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
791                  * (ypaprs(i,1)-ypplay(i,1))
792          END DO
793
794          ! Calculate the temperature et relative humidity at 2m and the wind at 10m
795          CALL stdlevvar(klon, knon, is_ter, zxli, &
796               yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
797               yts, yqsurf, yrugos, ypaprs(:,1), ypplay(:,1), &
798               yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
799         
800       END IF
801
802!****************************************************************************************
803!
804! 10) Switch selon current surface
805!     It is necessary to start with the continental surfaces because the ocean
806!     needs their run-off.
807!
808!****************************************************************************************
809       SELECT CASE(nsrf)
810     
811       CASE(is_ter)
812          ! ylwdown : to be removed, calculation is now done at land surface in surf_land
813          ylwdown(:)=0.0
814          DO i=1,knon
815             ylwdown(i)=lwdown_m(ni(i))
816          END DO
817          CALL surf_land(itap, dtime, date0, jour, knon, ni,&
818               rlon, rlat, &
819               debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
820               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
821               AcoefH, AcoefQ, BcoefH, BcoefQ, &
822               AcoefU, AcoefV, BcoefU, BcoefV, &
823               ypsref, yu1, yv1, yrugoro, pctsrf, &
824               ylwdown, yq2m, yt2m, &
825               ysnow, yqsol, yagesno, ytsoil, &
826               yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
827               yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
828               y_flux_u1, y_flux_v1 )
829               
830     
831       CASE(is_lic)
832          CALL surf_landice(itap, dtime, knon, ni, &
833               ysolsw, ysollw, yts, ypplay(:,1), &
834               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
835               AcoefH, AcoefQ, BcoefH, BcoefQ, &
836               AcoefU, AcoefV, BcoefU, BcoefV, &
837               ypsref, yu1, yv1, yrugoro, pctsrf, &
838               ysnow, yqsurf, yqsol, yagesno, &
839               ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
840               ytsurf_new, y_dflux_t, y_dflux_q, &
841               y_flux_u1, y_flux_v1)
842         
843       CASE(is_oce)
844          CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
845               yrugos, ywindsp, rmu0, yfder, yts, &
846               itap, dtime, jour, knon, ni, &
847               ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
848               AcoefH, AcoefQ, BcoefH, BcoefQ, &
849               AcoefU, AcoefV, BcoefU, BcoefV, &
850               ypsref, yu1, yv1, yrugoro, pctsrf, &
851               ysnow, yqsurf, yagesno, &
852               yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
853               ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
854               y_flux_u1, y_flux_v1)
855         
856       CASE(is_sic)
857          CALL surf_seaice( &
858               rlon, rlat, ysolsw, ysollw, yalb1, yfder, &
859               itap, dtime, jour, knon, ni, &
860               lafin, &
861               yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
862               AcoefH, AcoefQ, BcoefH, BcoefQ, &
863               AcoefU, AcoefV, BcoefU, BcoefV, &
864               ypsref, yu1, yv1, yrugoro, pctsrf, &
865               ysnow, yqsurf, yqsol, yagesno, ytsoil, &
866               yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
867               ytsurf_new, y_dflux_t, y_dflux_q, &
868               y_flux_u1, y_flux_v1)
869         
870
871       CASE DEFAULT
872          WRITE(lunout,*) 'Surface index = ', nsrf
873          abort_message = 'Surface index not valid'
874          CALL abort_gcm(modname,abort_message,1)
875       END SELECT
876
877
878!****************************************************************************************
879! 11) - Calcul the increment of surface temperature
880!
881!****************************************************************************************
882       y_d_ts(1:knon)   = ytsurf_new(1:knon) - yts(1:knon)
883 
884!****************************************************************************************
885!
886! 12) "La remontee" - "The uphill"
887!
888!  The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
889!  for X=H, Q, U and V, for all vertical levels.
890!
891!****************************************************************************************
892! H and Q
893       IF (ok_flux_surf) THEN
894          PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
895          y_flux_t1(:) =  fsens
896          y_flux_q1(:) =  flat/RLVTT
897          yfluxlat(:) =  flat
898       ELSE
899          y_flux_t1(:) =  yfluxsens(:)
900          y_flux_q1(:) = -yevap(:)
901       ENDIF
902
903       CALL climb_hq_up(knon, dtime, yt, yq, &
904            y_flux_q1, y_flux_t1, ypaprs, ypplay, &
905            y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))   
906       
907
908       CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
909            y_flux_u, y_flux_v, y_d_u, y_d_v)
910
911
912       DO j = 1, knon
913          y_dflux_t(j) = y_dflux_t(j) * ypct(j)
914          y_dflux_q(j) = y_dflux_q(j) * ypct(j)
915       ENDDO
916
917!****************************************************************************************
918! 13) Transform variables for output format :
919!     - Decompress
920!     - Multiply with pourcentage of current surface
921!     - Cumulate in global variable
922!
923!****************************************************************************************
924
925       tke(:,:,nsrf) = 0.
926       DO k = 1, klev
927          DO j = 1, knon
928             i = ni(j)
929             y_d_t(j,k)  = y_d_t(j,k) * ypct(j)
930             y_d_q(j,k)  = y_d_q(j,k) * ypct(j)
931             y_d_u(j,k)  = y_d_u(j,k) * ypct(j)
932             y_d_v(j,k)  = y_d_v(j,k) * ypct(j)
933
934             flux_t(i,k,nsrf) = y_flux_t(j,k)
935             flux_q(i,k,nsrf) = y_flux_q(j,k)
936             flux_u(i,k,nsrf) = y_flux_u(j,k)
937             flux_v(i,k,nsrf) = y_flux_v(j,k)
938
939             tke(i,k,nsrf)    = ytke(j,k)
940
941          ENDDO
942       ENDDO
943
944       evap(:,nsrf) = - flux_q(:,1,nsrf)
945       
946       alb1(:, nsrf) = 0.
947       alb2(:, nsrf) = 0.
948       snow(:, nsrf) = 0.
949       qsurf(:, nsrf) = 0.
950       rugos(:, nsrf) = 0.
951       fluxlat(:,nsrf) = 0.
952       DO j = 1, knon
953          i = ni(j)
954          d_ts(i,nsrf) = y_d_ts(j)
955          alb1(i,nsrf) = yalb1_new(j) 
956          alb2(i,nsrf) = yalb2_new(j)
957          snow(i,nsrf) = ysnow(j) 
958          qsurf(i,nsrf) = yqsurf(j)
959          rugos(i,nsrf) = yz0_new(j)
960          fluxlat(i,nsrf) = yfluxlat(j)
961          agesno(i,nsrf) = yagesno(j) 
962          cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
963          cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
964          dflux_t(i) = dflux_t(i) + y_dflux_t(j)
965          dflux_q(i) = dflux_q(i) + y_dflux_q(j)
966       END DO
967
968       DO k = 2, klev
969          DO j = 1, knon
970             i = ni(j)
971             zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)*ypct(j)
972          END DO
973       END DO
974
975       IF ( nsrf .EQ. is_ter ) THEN
976          DO j = 1, knon
977             i = ni(j)
978             qsol(i) = yqsol(j)
979          END DO
980       END IF
981       
982       ftsoil(:,:,nsrf) = 0.
983       DO k = 1, nsoilmx
984          DO j = 1, knon
985             i = ni(j)
986             ftsoil(i, k, nsrf) = ytsoil(j,k)
987          END DO
988       END DO
989       
990       
991       DO k = 1, klev
992          DO j = 1, knon
993             i = ni(j)
994             d_t(i,k) = d_t(i,k) + y_d_t(j,k)
995             d_q(i,k) = d_q(i,k) + y_d_q(j,k)
996             d_u(i,k) = d_u(i,k) + y_d_u(j,k)
997             d_v(i,k) = d_v(i,k) + y_d_v(j,k)
998          END DO
999       END DO
1000
1001!****************************************************************************************
1002! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
1003!     Call HBTM
1004!
1005!****************************************************************************************
1006       t2m(:,nsrf)    = 0.
1007       q2m(:,nsrf)    = 0.
1008       u10m(:,nsrf)   = 0.
1009       v10m(:,nsrf)   = 0.
1010
1011       pblh(:,nsrf)   = 0.        ! Hauteur de couche limite
1012       plcl(:,nsrf)   = 0.        ! Niveau de condensation de la CLA
1013       capCL(:,nsrf)  = 0.        ! CAPE de couche limite
1014       oliqCL(:,nsrf) = 0.        ! eau_liqu integree de couche limite
1015       cteiCL(:,nsrf) = 0.        ! cloud top instab. crit. couche limite
1016       pblt(:,nsrf)   = 0.        ! T a la Hauteur de couche limite
1017       therm(:,nsrf)  = 0.
1018       trmb1(:,nsrf)  = 0.        ! deep_cape
1019       trmb2(:,nsrf)  = 0.        ! inhibition
1020       trmb3(:,nsrf)  = 0.        ! Point Omega
1021
1022#undef T2m     
1023#define T2m     
1024#ifdef T2m
1025! Calculations of diagnostic t,q at 2m and u, v at 10m
1026
1027       DO j=1, knon
1028          i = ni(j)
1029          uzon(j) = yu(j,1) + y_d_u(j,1)
1030          vmer(j) = yv(j,1) + y_d_v(j,1)
1031          tair1(j) = yt(j,1) + y_d_t(j,1)
1032          qair1(j) = yq(j,1) + y_d_q(j,1)
1033          zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
1034               * (ypaprs(j,1)-ypplay(j,1))
1035          tairsol(j) = yts(j) + y_d_ts(j)
1036          rugo1(j) = yrugos(j)
1037          IF(nsrf.EQ.is_oce) THEN
1038             rugo1(j) = rugos(i,nsrf)
1039          ENDIF
1040          psfce(j)=ypaprs(j,1)
1041          patm(j)=ypplay(j,1)
1042          qairsol(j) = yqsurf(j)
1043       END DO
1044       
1045
1046! Calculate the temperature et relative humidity at 2m and the wind at 10m
1047       CALL stdlevvar(klon, knon, nsrf, zxli, &
1048            uzon, vmer, tair1, qair1, zgeo1, &
1049            tairsol, qairsol, rugo1, psfce, patm, &
1050            yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
1051
1052       DO j=1, knon
1053          i = ni(j)
1054          t2m(i,nsrf)=yt2m(j)
1055          q2m(i,nsrf)=yq2m(j)
1056         
1057          ! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
1058          u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)**2+vmer(j)**2)
1059          v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)**2+vmer(j)**2)
1060       END DO
1061
1062!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
1063!IM Ajoute dependance type surface
1064       IF (thermcep) THEN
1065          DO j = 1, knon
1066             i=ni(j)
1067             zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
1068             zx_qs1  = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
1069             zx_qs1  = MIN(0.5,zx_qs1)
1070             zcor1   = 1./(1.-RETV*zx_qs1)
1071             zx_qs1  = zx_qs1*zcor1
1072             
1073             rh2m(i)   = rh2m(i)   + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
1074             qsat2m(i) = qsat2m(i) + zx_qs1  * pctsrf(i,nsrf)
1075          END DO
1076       END IF
1077
1078       CALL HBTM(knon, ypaprs, ypplay, &
1079            yt2m,yt10m,yq2m,yq10m,yustar, &
1080            y_flux_t,y_flux_q,yu,yv,yt,yq, &
1081            ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
1082            ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
1083       
1084       DO j=1, knon
1085          i = ni(j)
1086          pblh(i,nsrf)   = ypblh(j)
1087          plcl(i,nsrf)   = ylcl(j)
1088          capCL(i,nsrf)  = ycapCL(j)
1089          oliqCL(i,nsrf) = yoliqCL(j)
1090          cteiCL(i,nsrf) = ycteiCL(j)
1091          pblT(i,nsrf)   = ypblT(j)
1092          therm(i,nsrf)  = ytherm(j)
1093          trmb1(i,nsrf)  = ytrmb1(j)
1094          trmb2(i,nsrf)  = ytrmb2(j)
1095          trmb3(i,nsrf)  = ytrmb3(j)
1096       END DO
1097       
1098#else
1099! T2m not defined
1100! No calculation
1101       PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
1102#endif
1103
1104!****************************************************************************************
1105! 15) End of loop over different surfaces
1106!
1107!****************************************************************************************
1108    END DO loop_nbsrf
1109
1110!****************************************************************************************
1111! 16) Calculate the mean value over all sub-surfaces for som variables
1112!
1113!****************************************************************************************
1114   
1115    zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
1116    zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
1117    DO nsrf = 1, nbsrf
1118       DO k = 1, klev
1119          DO i = 1, klon
1120             zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
1121             zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
1122             zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
1123             zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
1124          END DO
1125       END DO
1126    END DO
1127
1128    DO i = 1, klon
1129       zxsens(i)     = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1130       zxevap(i)     = - zxfluxq(i,1) ! flux d'evaporation au sol
1131       fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
1132    ENDDO
1133   
1134!
1135! Incrementer la temperature du sol
1136!
1137    zxtsol(:) = 0.0  ; zxfluxlat(:) = 0.0
1138    zt2m(:) = 0.0    ; zq2m(:) = 0.0
1139    zu10m(:) = 0.0   ; zv10m(:) = 0.0
1140    s_pblh(:) = 0.0  ; s_plcl(:) = 0.0
1141    s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
1142    s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
1143    s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
1144    s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
1145   
1146   
1147    DO nsrf = 1, nbsrf
1148       DO i = 1, klon         
1149          ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
1150         
1151          wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
1152               + flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1153          wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
1154               pctsrf(i,nsrf)
1155
1156          zxtsol(i)    = zxtsol(i)    + ts(i,nsrf)      * pctsrf(i,nsrf)
1157          zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
1158         
1159          zt2m(i)  = zt2m(i)  + t2m(i,nsrf)  * pctsrf(i,nsrf)
1160          zq2m(i)  = zq2m(i)  + q2m(i,nsrf)  * pctsrf(i,nsrf)
1161          zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
1162          zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
1163
1164          s_pblh(i)   = s_pblh(i)   + pblh(i,nsrf)  * pctsrf(i,nsrf)
1165          s_plcl(i)   = s_plcl(i)   + plcl(i,nsrf)  * pctsrf(i,nsrf)
1166          s_capCL(i)  = s_capCL(i)  + capCL(i,nsrf) * pctsrf(i,nsrf)
1167          s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
1168          s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
1169          s_pblT(i)   = s_pblT(i)   + pblT(i,nsrf)  * pctsrf(i,nsrf)
1170          s_therm(i)  = s_therm(i)  + therm(i,nsrf) * pctsrf(i,nsrf)
1171          s_trmb1(i)  = s_trmb1(i)  + trmb1(i,nsrf) * pctsrf(i,nsrf)
1172          s_trmb2(i)  = s_trmb2(i)  + trmb2(i,nsrf) * pctsrf(i,nsrf)
1173          s_trmb3(i)  = s_trmb3(i)  + trmb3(i,nsrf) * pctsrf(i,nsrf)
1174       END DO
1175    END DO
1176
1177    IF (check) THEN
1178       amn=MIN(ts(1,is_ter),1000.)
1179       amx=MAX(ts(1,is_ter),-1000.)
1180       DO i=2, klon
1181          amn=MIN(ts(i,is_ter),amn)
1182          amx=MAX(ts(i,is_ter),amx)
1183       ENDDO
1184       PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
1185    ENDIF
1186
1187!jg ?
1188!!$!
1189!!$! If a sub-surface does not exsist for a grid point, the mean value for all
1190!!$! sub-surfaces is distributed.
1191!!$!
1192!!$    DO nsrf = 1, nbsrf
1193!!$       DO i = 1, klon
1194!!$          IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
1195!!$             ts(i,nsrf)     = zxtsol(i)
1196!!$             t2m(i,nsrf)    = zt2m(i)
1197!!$             q2m(i,nsrf)    = zq2m(i)
1198!!$             u10m(i,nsrf)   = zu10m(i)
1199!!$             v10m(i,nsrf)   = zv10m(i)
1200!!$
1201!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
1202!!$             pblh(i,nsrf)   = s_pblh(i)
1203!!$             plcl(i,nsrf)   = s_plcl(i)
1204!!$             capCL(i,nsrf)  = s_capCL(i)
1205!!$             oliqCL(i,nsrf) = s_oliqCL(i)
1206!!$             cteiCL(i,nsrf) = s_cteiCL(i)
1207!!$             pblT(i,nsrf)   = s_pblT(i)
1208!!$             therm(i,nsrf)  = s_therm(i)
1209!!$             trmb1(i,nsrf)  = s_trmb1(i)
1210!!$             trmb2(i,nsrf)  = s_trmb2(i)
1211!!$             trmb3(i,nsrf)  = s_trmb3(i)
1212!!$          ENDIF
1213!!$       ENDDO
1214!!$    ENDDO
1215
1216
1217    DO i = 1, klon
1218       fder(i) = - 4.0*RSIGMA*zxtsol(i)**3
1219    ENDDO
1220   
1221    zxqsurf(:) = 0.0
1222    zxsnow(:)  = 0.0
1223    DO nsrf = 1, nbsrf
1224       DO i = 1, klon
1225          zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
1226          zxsnow(i)  = zxsnow(i)  + snow(i,nsrf)  * pctsrf(i,nsrf)
1227       END DO
1228    END DO
1229
1230! Premier niveau de vent sortie dans physiq.F
1231    zu1(:) = u(:,1)
1232    zv1(:) = v(:,1)
1233
1234! Some of the module declared variables are returned for printing in physiq.F
1235    qsol_d(:)     = qsol(:)
1236    evap_d(:,:)   = evap(:,:)
1237    rugos_d(:,:)  = rugos(:,:)
1238    agesno_d(:,:) = agesno(:,:)
1239
1240
1241  END SUBROUTINE pbl_surface
1242!
1243!****************************************************************************************
1244!
1245  SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
1246       evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
1247
1248    INCLUDE "indicesol.h"
1249    INCLUDE "dimsoil.h"
1250
1251! Ouput variables
1252!****************************************************************************************
1253    REAL, DIMENSION(klon), INTENT(OUT)                 :: qsol_rst
1254    REAL, DIMENSION(klon), INTENT(OUT)                 :: fder_rst
1255    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: snow_rst
1256    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: qsurf_rst
1257    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: evap_rst
1258    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: rugos_rst
1259    REAL, DIMENSION(klon, nbsrf), INTENT(OUT)          :: agesno_rst
1260    REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
1261
1262 
1263!****************************************************************************************
1264! Return module variables for writing to restart file
1265!
1266!****************************************************************************************   
1267    qsol_rst(:)       = qsol(:)
1268    fder_rst(:)       = fder(:)
1269    snow_rst(:,:)     = snow(:,:)
1270    qsurf_rst(:,:)    = qsurf(:,:)
1271    evap_rst(:,:)     = evap(:,:)
1272    rugos_rst(:,:)    = rugos(:,:)
1273    agesno_rst(:,:)   = agesno(:,:)
1274    ftsoil_rst(:,:,:) = ftsoil(:,:,:)
1275
1276!****************************************************************************************
1277! Deallocate module variables
1278!
1279!****************************************************************************************
1280    DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
1281
1282  END SUBROUTINE pbl_surface_final
1283
1284!****************************************************************************************
1285!
1286  SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, tke)
1287
1288    ! Give default values where new fraction has appread
1289
1290    INCLUDE "indicesol.h"
1291    INCLUDE "dimsoil.h"
1292    INCLUDE "clesphys.h"
1293    INCLUDE "compbl.h"
1294
1295! Input variables
1296!****************************************************************************************
1297    INTEGER, INTENT(IN)                     :: itime
1298    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
1299
1300! InOutput variables
1301!****************************************************************************************
1302    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: tsurf
1303    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: alb1, alb2
1304    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT)        :: u10m, v10m
1305    REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
1306
1307! Local variables
1308!****************************************************************************************
1309    INTEGER           :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
1310    CHARACTER(len=80) :: abort_message
1311    CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
1312    INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
1313!
1314! All at once !!
1315!****************************************************************************************
1316   
1317    DO nsrf = 1, nbsrf
1318       ! First decide complement sub-surfaces
1319       SELECT CASE (nsrf)
1320       CASE(is_oce)
1321          nsrf_comp1=is_sic
1322          nsrf_comp2=is_ter
1323          nsrf_comp3=is_lic
1324       CASE(is_sic)
1325          nsrf_comp1=is_oce
1326          nsrf_comp2=is_ter
1327          nsrf_comp3=is_lic
1328       CASE(is_ter)
1329          nsrf_comp1=is_lic
1330          nsrf_comp2=is_oce
1331          nsrf_comp3=is_sic
1332       CASE(is_lic)
1333          nsrf_comp1=is_ter
1334          nsrf_comp2=is_oce
1335          nsrf_comp3=is_sic
1336       END SELECT
1337
1338       ! Initialize all new fractions
1339       DO i=1, klon
1340          IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
1341             
1342             IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
1343                ! Use the complement sub-surface, keeping the continents unchanged
1344                qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
1345                evap(i,nsrf)  = evap(i,nsrf_comp1)
1346                rugos(i,nsrf) = rugos(i,nsrf_comp1)
1347                tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
1348                alb1(i,nsrf)  = alb1(i,nsrf_comp1)
1349                alb2(i,nsrf)  = alb2(i,nsrf_comp1)
1350                u10m(i,nsrf)  = u10m(i,nsrf_comp1)
1351                v10m(i,nsrf)  = v10m(i,nsrf_comp1)
1352                if (iflag_pbl > 1) then
1353                 tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
1354                endif
1355                mfois(nsrf) = mfois(nsrf) + 1
1356             ELSE
1357                ! The continents have changed. The new fraction receives the mean sum of the existent fractions
1358                qsurf(i,nsrf) = qsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
1359                evap(i,nsrf)  = evap(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1360                rugos(i,nsrf) = rugos(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
1361                tsurf(i,nsrf) = tsurf(i,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
1362                alb1(i,nsrf)  = alb1(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1363                alb2(i,nsrf)  = alb2(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1364                u10m(i,nsrf)  = u10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1365                v10m(i,nsrf)  = v10m(i,nsrf_comp2) *pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) *pctsrf_old(i,nsrf_comp3)
1366                if (iflag_pbl > 1) then
1367                 tke(i,:,nsrf) = tke(i,:,nsrf_comp2)*pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)*pctsrf_old(i,nsrf_comp3)
1368                endif
1369           
1370                ! Security abort. This option has never been tested. To test, comment the following line.
1371!                abort_message='The fraction of the continents have changed!'
1372!                CALL abort_gcm(modname,abort_message,1)
1373                nfois(nsrf) = nfois(nsrf) + 1
1374             END IF
1375             snow(i,nsrf)     = 0.
1376             agesno(i,nsrf)   = 0.
1377             ftsoil(i,:,nsrf) = tsurf(i,nsrf)
1378          ELSE
1379             pfois(nsrf) = pfois(nsrf)+ 1
1380          END IF
1381       END DO
1382       
1383    END DO
1384
1385  END SUBROUTINE pbl_surface_newfrac
1386
1387
1388!****************************************************************************************
1389
1390
1391END MODULE pbl_surface_mod
Note: See TracBrowser for help on using the repository browser.