Context Navigation

pbl_surface_mod.F90 @ 2031

Last change on this file since 2031 was 1932, checked in by lguez, 11 years ago

Declaration of variables appearing in array bounds must be specified
before the declaration of the array.

Variables used in FCTTRE.h are declared in YOMCST.h and YOETHF.h so
YOMCST.h and YOETHF.h must be included before FCTTRE.h.

(There were compilation errors with gfortran and debugging options.)

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

File size: 64.6 KB

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

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

Download in other formats:

Original Format