Context Navigation

pbl_surface_mod.F90 @ 5447

Last change on this file since 5447 was 1539, checked in by musat, 14 years ago

Ajouts CFMIP2/CMIP5

6eme fichier de sortie "stations" histstn.nc qui necessite 2 fichiers (voir DefLists?): npCFMIP_param.data(_*) contenant le nombre de points (120 pour simulations AMIP, 73 pour aqua) pointlocations.txt(_*) contenant 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

Ajout flag LOGICAL lCOSP necessaire pour sortir un fichier stations
IM

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

File size: 58.8 KB

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

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

Download in other formats:

Original Format