Context Navigation

pbl_surface_mod.F90 @ 1563

Last change on this file since 1563 was 1555, checked in by jghattas, 13 years ago
Ajout du calcul de la temperature sol pour le 1D si flux forces(Nicolas Rochetin).
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 59.1 KB

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

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

Download in other formats:

Original Format