Context Navigation

pbl_surface_mod.F90 @ 2370

Last change on this file since 2370 was 1534, checked in by musat, 13 years ago

Ajouts CFMIP2/CMIP5

6eme fichier de sortie "stations" histstn.nc qui necessite 2 fichiers: PARAM/npCFMIP_param.data contenant le nombre de points (120 pour simulations AMIP, 73 pour aqua) PARAM/pointlocations.txt contenat le numero, les coordonnees (lon,lat) et le nom de chaque station
flag LOGICAL dans tous les appels histwrite_phy pour pouvoir sortir le fichier histstn.nc

NB: 1) les flags de type phys_ que l'on met dans le physiq.def_L* pour ajouter plus de sorties

necessitent dorenavant 6 valeurs, la 6eme correspondant au fichier histstn.nc

2) par defaut le fichier histstn.nc ne sort pas; pour le sortir ajouter les lignes suivantes

dans physiq.def_L*

### Type de fichier : global (n) ou stations (y)
phys_out_filestations = n n n n n y

introduction de 2 jeux de flags pour les taux des GES; taux actuels avec suffixes _act, taux futurs avec "_per" avec 2 appels au rayonnement si taux "_per" different des taux "_act" (utiles pour diags. CFMIP 4CO2)
flags "betaCRF" pour calculs CRF pour experiences sensibilite proprietes optiques eau liquide nuageuse avec initialisations par defaut; sinon besoin de fichier beta_crf.data

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 58.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format