Context Navigation

pbl_surface_mod.F90 @ 3809

Last change on this file since 3809 was 1942, checked in by jghattas, 11 years ago
Passing the cosine of solar zenith angle from pbl_surface into orchidee. Matthew McGrath?(LSCE)
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 61.5 KB

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

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

Download in other formats:

Original Format