Context Navigation

pbl_surface_mod.F90 @ 2153

Last change on this file since 2153 was 2126, checked in by fhourdin, 10 years ago

Introduction du cas Dice couplé atmosphere/surface
+ nouveau paramètre de contrôle f_ri_cd_min, seuil minimum
sur la fonction f(Ri) en facteur du coefficient de traîné neutre.
Par défaut : f_ri_cd_min=0.1 (comme avant)

Introduction of the coupled atmosphere/surface dice case.
+ new parameter f_ri_cd_min, minimum threshold on the f(Ri) fonction
that multiply the neutral drag coefficient.

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: 65.8 KB

Rev	Line
[781]	1	!
[1279]	2	! $Id: pbl_surface_mod.F90 2126 2014-10-15 00:03:57Z lguez $
	3	!
[781]	4	MODULE pbl_surface_mod
	5	!
	6	! Planetary Boundary Layer and Surface module
	7	!
	8	! This module manage the calculation of turbulent diffusion in the boundary layer
	9	! and all interactions towards the differents sub-surfaces.
	10	!
	11	!
	12	USE dimphy
	13	USE mod_phys_lmdz_para, ONLY : mpi_size
	14	USE ioipsl
[996]	15	USE surface_data, ONLY : type_ocean, ok_veget
[781]	16	USE surf_land_mod, ONLY : surf_land
	17	USE surf_landice_mod, ONLY : surf_landice
	18	USE surf_ocean_mod, ONLY : surf_ocean
	19	USE surf_seaice_mod, ONLY : surf_seaice
	20	USE cpl_mod, ONLY : gath2cpl
	21	USE climb_hq_mod, ONLY : climb_hq_down, climb_hq_up
	22	USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up
	23	USE coef_diff_turb_mod, ONLY : coef_diff_turb
[1403]	24	USE control_mod
[781]	25
[1403]	26
[781]	27	IMPLICIT NONE
	28
	29	! Declaration of variables saved in restart file
[888]	30	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: qsol ! water height in the soil (mm)
[781]	31	!$OMP THREADPRIVATE(qsol)
[888]	32	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: fder ! flux drift
[781]	33	!$OMP THREADPRIVATE(fder)
[888]	34	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: snow ! snow at surface
[781]	35	!$OMP THREADPRIVATE(snow)
[888]	36	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qsurf ! humidity at surface
[781]	37	!$OMP THREADPRIVATE(qsurf)
[888]	38	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: evap ! evaporation at surface
[781]	39	!$OMP THREADPRIVATE(evap)
[888]	40	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: rugos ! rugosity at surface (m)
[781]	41	!$OMP THREADPRIVATE(rugos)
[888]	42	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: agesno ! age of snow at surface
[781]	43	!$OMP THREADPRIVATE(agesno)
[1787]	44	! Correction pour le cas AMMA (PRIVATE)
	45	REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ftsoil ! soil temperature
[781]	46	!$OMP THREADPRIVATE(ftsoil)
	47
	48	CONTAINS
	49	!
	50	!****************************************************************************************
	51	!
	52	SUBROUTINE pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst,&
	53	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
	54
	55	! This routine should be called after the restart file has been read.
	56	! This routine initialize the restart variables and does some validation tests
	57	! for the index of the different surfaces and tests the choice of type of ocean.
	58
[1785]	59	USE indice_sol_mod
	60
[781]	61	INCLUDE "dimsoil.h"
	62	INCLUDE "iniprint.h"
	63
	64	! Input variables
	65	!****************************************************************************************
	66	REAL, DIMENSION(klon), INTENT(IN) :: qsol_rst
	67	REAL, DIMENSION(klon), INTENT(IN) :: fder_rst
	68	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: snow_rst
	69	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: qsurf_rst
	70	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: evap_rst
	71	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: rugos_rst
	72	REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: agesno_rst
	73	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(IN) :: ftsoil_rst
	74
	75
	76	! Local variables
	77	!****************************************************************************************
	78	INTEGER :: ierr
	79	CHARACTER(len=80) :: abort_message
	80	CHARACTER(len = 20) :: modname = 'pbl_surface_init'
	81
	82
	83	!****************************************************************************************
	84	! Allocate and initialize module variables with fields read from restart file.
	85	!
	86	!****************************************************************************************
	87	ALLOCATE(qsol(klon), stat=ierr)
	88	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
	89
	90	ALLOCATE(fder(klon), stat=ierr)
	91	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
	92
	93	ALLOCATE(snow(klon,nbsrf), stat=ierr)
	94	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
	95
	96	ALLOCATE(qsurf(klon,nbsrf), stat=ierr)
	97	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
	98
	99	ALLOCATE(evap(klon,nbsrf), stat=ierr)
	100	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
	101
	102	ALLOCATE(rugos(klon,nbsrf), stat=ierr)
	103	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
	104
	105	ALLOCATE(agesno(klon,nbsrf), stat=ierr)
	106	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
	107
	108	ALLOCATE(ftsoil(klon,nsoilmx,nbsrf), stat=ierr)
	109	IF (ierr /= 0) CALL abort_gcm('pbl_surface_init', 'pb in allocation',1)
	110
	111
	112	qsol(:) = qsol_rst(:)
	113	fder(:) = fder_rst(:)
	114	snow(:,:) = snow_rst(:,:)
	115	qsurf(:,:) = qsurf_rst(:,:)
	116	evap(:,:) = evap_rst(:,:)
	117	rugos(:,:) = rugos_rst(:,:)
	118	agesno(:,:) = agesno_rst(:,:)
	119	ftsoil(:,:,:) = ftsoil_rst(:,:,:)
	120
	121
	122	!****************************************************************************************
	123	! Test for sub-surface indices
	124	!
	125	!****************************************************************************************
	126	IF (is_ter /= 1) THEN
	127	WRITE(lunout,)" Warning *"
	128	WRITE(lunout,*)" is_ter n'est pas le premier surface, is_ter = ",is_ter
	129	WRITE(lunout,*)"or on doit commencer par les surfaces continentales"
	130	abort_message="voir ci-dessus"
	131	CALL abort_gcm(modname,abort_message,1)
	132	ENDIF
	133
	134	IF ( is_oce > is_sic ) THEN
	135	WRITE(lunout,)' Warning *'
	136	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
	137	WRITE(lunout,*)' l''ocean doit etre traite avant la banquise'
	138	WRITE(lunout,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
	139	abort_message='voir ci-dessus'
	140	CALL abort_gcm(modname,abort_message,1)
	141	ENDIF
	142
	143	IF ( is_lic > is_sic ) THEN
	144	WRITE(lunout,)' Warning *'
	145	WRITE(lunout,*)' Pour des raisons de sequencement dans le code'
	146	WRITE(lunout,*)' la glace contineltalle doit etre traite avant la glace de mer'
	147	WRITE(lunout,*)' or is_lic = ',is_lic, '> is_sic = ',is_sic
	148	abort_message='voir ci-dessus'
	149	CALL abort_gcm(modname,abort_message,1)
	150	ENDIF
	151
	152	!****************************************************************************************
	153	! Validation of ocean mode
	154	!
	155	!****************************************************************************************
	156
[996]	157	IF (type_ocean /= 'slab ' .AND. type_ocean /= 'force ' .AND. type_ocean /= 'couple') THEN
[1064]	158	WRITE(lunout,)' Warning *'
	159	WRITE(lunout,*)'Option couplage pour l''ocean = ', type_ocean
	160	abort_message='option pour l''ocean non valable'
	161	CALL abort_gcm(modname,abort_message,1)
[781]	162	ENDIF
	163
	164	END SUBROUTINE pbl_surface_init
	165	!
	166	!****************************************************************************************
	167	!
	168
	169	SUBROUTINE pbl_surface( &
	170	dtime, date0, itap, jour, &
	171	debut, lafin, &
	172	rlon, rlat, rugoro, rmu0, &
[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
[2126]	321	REAL, DIMENSION(klon), INTENT(OUT) :: alb3_lic
[1865]	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
[2126]	477	REAL :: vent
[781]	478
	479	! For debugging with IOIPSL
	480	INTEGER, DIMENSION(iim*(jjm+1)) :: ndexbg
	481	REAL :: zjulian
	482	REAL, DIMENSION(klon) :: tabindx
	483	REAL, DIMENSION(iim,jjm+1) :: zx_lon, zx_lat
	484	REAL, DIMENSION(iim,jjm+1) :: debugtab
	485
	486
[888]	487	REAL, DIMENSION(klon,nbsrf) :: pblh ! height of the planetary boundary layer
	488	REAL, DIMENSION(klon,nbsrf) :: plcl ! condensation level
[781]	489	REAL, DIMENSION(klon,nbsrf) :: capCL
	490	REAL, DIMENSION(klon,nbsrf) :: oliqCL
	491	REAL, DIMENSION(klon,nbsrf) :: cteiCL
	492	REAL, DIMENSION(klon,nbsrf) :: pblT
	493	REAL, DIMENSION(klon,nbsrf) :: therm
[888]	494	REAL, DIMENSION(klon,nbsrf) :: trmb1 ! deep cape
	495	REAL, DIMENSION(klon,nbsrf) :: trmb2 ! inhibition
	496	REAL, DIMENSION(klon,nbsrf) :: trmb3 ! point Omega
[1067]	497	REAL, DIMENSION(klon,nbsrf) :: zx_rh2m, zx_qsat2m
[996]	498	REAL, DIMENSION(klon,nbsrf) :: zx_t1
[888]	499	REAL, DIMENSION(klon, nbsrf) :: alb ! mean albedo for whole SW interval
	500	REAL, DIMENSION(klon) :: ylwdown ! jg : temporary (ysollwdown)
[781]	501
[996]	502	REAL :: zx_qs1, zcor1, zdelta1
[781]	503
[1865]	504	! Martin
	505	REAL, DIMENSION(klon, nbsrf) :: sollwd ! net longwave radiation at surface
	506	REAL, DIMENSION(klon) :: ysollwd
	507	REAL, DIMENSION(klon) :: ytoice
	508	REAL, DIMENSION(klon) :: ysnowhgt, yqsnow, ysissnow, yrunoff
	509	REAL, DIMENSION(klon) :: yzsig
	510	REAL, DIMENSION(klon,klev) :: ypphi
	511	REAL, DIMENSION(klon) :: ycldt
	512	REAL, DIMENSION(klon) :: yrmu0
	513	! Martin
	514
[781]	515	!****************************************************************************************
[1894]	516
[781]	517	! End of declarations
	518	!****************************************************************************************
	519
	520
	521	!****************************************************************************************
	522	! 1) Initialisation and validation tests
	523	! Only done first time entering this subroutine
	524	!
	525	!****************************************************************************************
	526
	527	IF (first_call) THEN
	528	first_call=.FALSE.
	529
[1282]	530	! Initialize ok_flux_surf (for 1D model)
[1403]	531	if (klon>1) ok_flux_surf=.FALSE.
[1282]	532
[781]	533	! Initilize debug IO
	534	IF (debugindex .AND. mpi_size==1) THEN
	535	! initialize IOIPSL output
	536	idayref = day_ini
	537	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
	538	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
	539	DO i = 1, iim
	540	zx_lon(i,1) = rlon(i+1)
	541	zx_lon(i,jjm+1) = rlon(i+1)
	542	ENDDO
	543	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
	544	CALL histbeg("sous_index", iim,zx_lon(:,1),jjm+1,zx_lat(1,:), &
	545	1,iim,1,jjm+1, &
	546	itau_phy,zjulian,dtime,nhoridbg,nidbg)
	547	! no vertical axis
	548	cl_surf(1)='ter'
	549	cl_surf(2)='lic'
	550	cl_surf(3)='oce'
	551	cl_surf(4)='sic'
	552	DO nsrf=1,nbsrf
	553	CALL histdef(nidbg, cl_surf(nsrf),cl_surf(nsrf), "-",iim, &
	554	jjm+1,nhoridbg, 1, 1, 1, -99, 32, "inst", dtime,dtime)
	555	END DO
	556
	557	CALL histend(nidbg)
	558	CALL histsync(nidbg)
	559
	560	END IF
	561
	562	ENDIF
	563
	564	!****************************************************************************************
[889]	565	! Force soil water content to qsol0 if qsol0>0 and VEGET=F (use bucket
	566	! instead of ORCHIDEE)
[1894]	567	IF (qsol0>=0.) THEN
[1067]	568	PRINT*,'WARNING : On impose qsol=',qsol0
[889]	569	qsol(:)=qsol0
[1067]	570	ENDIF
[889]	571	!****************************************************************************************
	572
	573	!****************************************************************************************
[781]	574	! 2) Initialization to zero
	575	! Done for all local variables that will be compressed later
	576	! and argument with INTENT(OUT)
	577	!****************************************************************************************
	578	cdragh = 0.0 ; cdragm = 0.0 ; dflux_t = 0.0 ; dflux_q = 0.0
[1067]	579	ypct = 0.0 ; yts = 0.0 ; ysnow = 0.0
[888]	580	zv1 = 0.0 ; yqsurf = 0.0 ; yalb1 = 0.0 ; yalb2 = 0.0
[781]	581	yrain_f = 0.0 ; ysnow_f = 0.0 ; yfder = 0.0 ; ysolsw = 0.0
[888]	582	ysollw = 0.0 ; yrugos = 0.0 ; yu1 = 0.0
	583	yv1 = 0.0 ; ypaprs = 0.0 ; ypplay = 0.0
[781]	584	ydelp = 0.0 ; yu = 0.0 ; yv = 0.0 ; yt = 0.0
[996]	585	yq = 0.0 ; y_dflux_t = 0.0 ; y_dflux_q = 0.0
[1067]	586	yrugoro = 0.0 ; ywindsp = 0.0
[781]	587	d_ts = 0.0 ; yfluxlat=0.0 ; flux_t = 0.0 ; flux_q = 0.0
	588	flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0
[1761]	589	d_t_diss= 0.0 ;d_u = 0.0 ; d_v = 0.0 ; yqsol = 0.0
[878]	590	ytherm = 0.0 ; ytke=0.
[1865]	591	! Martin
	592	ysnowhgt = 0.0; yqsnow = 0.0 ; yrunoff = 0.0 ; ytoice =0.0
	593	yalb3_new = 0.0 ; ysissnow = 0.0 ; ysollwd = 0.0
	594	ypphi = 0.0 ; ycldt = 0.0 ; yrmu0 = 0.0
	595	! Martin
[1067]	596
[1761]	597	tke(:,:,is_ave)=0.
	598	IF (iflag_pbl<20.or.iflag_pbl>=30) THEN
	599	zcoefh(:,:,:) = 0.0
	600	zcoefh(:,1,:) = 999999. ! zcoefh(:,k=1) should never be used
	601	zcoefm(:,:,:) = 0.0
	602	zcoefm(:,1,:) = 999999. !
	603	ELSE
	604	zcoefm(:,:,is_ave)=0.
	605	zcoefh(:,:,is_ave)=0.
	606	ENDIF
[781]	607	ytsoil = 999999.
	608
[1064]	609	rh2m(:) = 0.
	610	qsat2m(:) = 0.
[781]	611	!****************************************************************************************
	612	! 3) - Calculate pressure thickness of each layer
	613	! - Calculate the wind at first layer
[888]	614	! - Mean calculations of albedo
	615	! - Calculate net radiance at sub-surface
[781]	616	!****************************************************************************************
	617	DO k = 1, klev
	618	DO i = 1, klon
	619	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	620	ENDDO
	621	ENDDO
	622
	623	!****************************************************************************************
	624	! Test for rugos........ from physiq.. A la fin plutot???
[888]	625	!
[781]	626	!****************************************************************************************
	627
	628	zxrugs(:) = 0.0
	629	DO nsrf = 1, nbsrf
	630	DO i = 1, klon
	631	rugos(i,nsrf) = MAX(rugos(i,nsrf),0.000015)
	632	zxrugs(i) = zxrugs(i) + rugos(i,nsrf)*pctsrf(i,nsrf)
	633	ENDDO
	634	ENDDO
	635
[888]	636	! Mean calculations of albedo
	637	!
	638	! Albedo at sub-surface
	639	! * alb1 : albedo in visible SW interval
	640	! * alb2 : albedo in near infrared SW interval
	641	! * alb : mean albedo for whole SW interval
	642	!
	643	! Mean albedo for grid point
	644	! * alb1_m : albedo in visible SW interval
	645	! * alb2_m : albedo in near infrared SW interval
	646	! * alb_m : mean albedo at whole SW interval
	647
	648	alb1_m(:) = 0.0
	649	alb2_m(:) = 0.0
[781]	650	DO nsrf = 1, nbsrf
	651	DO i = 1, klon
[888]	652	alb1_m(i) = alb1_m(i) + alb1(i,nsrf) * pctsrf(i,nsrf)
	653	alb2_m(i) = alb2_m(i) + alb2(i,nsrf) * pctsrf(i,nsrf)
[781]	654	ENDDO
	655	ENDDO
	656
[888]	657	! We here suppose the fraction f1 of incoming radiance of visible radiance
	658	! as a fraction of all shortwave radiance
[1069]	659	f1 = 0.5
	660	! f1 = 1 ! put f1=1 to recreate old calculations
[781]	661
[888]	662	DO nsrf = 1, nbsrf
	663	DO i = 1, klon
	664	alb(i,nsrf) = f1alb1(i,nsrf) + (1-f1)alb2(i,nsrf)
	665	ENDDO
	666	ENDDO
[781]	667
[888]	668	DO i = 1, klon
	669	alb_m(i) = f1alb1_m(i) + (1-f1)alb2_m(i)
	670	END DO
	671
	672	! Calculation of mean temperature at surface grid points
[781]	673	ztsol(:) = 0.0
	674	DO nsrf = 1, nbsrf
	675	DO i = 1, klon
	676	ztsol(i) = ztsol(i) + ts(i,nsrf)*pctsrf(i,nsrf)
	677	ENDDO
	678	ENDDO
	679
[888]	680	! Linear distrubution on sub-surface of long- and shortwave net radiance
[781]	681	DO nsrf = 1, nbsrf
	682	DO i = 1, klon
	683	sollw(i,nsrf) = sollw_m(i) + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ts(i,nsrf))
[1865]	684	! Martin
	685	sollwd(i,nsrf)= sollwd_m(i)
	686	! Martin
[888]	687	solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
[781]	688	ENDDO
	689	ENDDO
	690
	691
[888]	692	! Downwelling longwave radiation at mean surface
	693	lwdown_m(:) = 0.0
[781]	694	DO i = 1, klon
[888]	695	lwdown_m(i) = sollw_m(i) + RSIGMAztsol(i)*4
[781]	696	ENDDO
	697
	698	!****************************************************************************************
	699	! 4) Loop over different surfaces
	700	!
[996]	701	! Only points containing a fraction of the sub surface will be threated.
[781]	702	!
	703	!****************************************************************************************
[1064]	704
[781]	705	loop_nbsrf: DO nsrf = 1, nbsrf
	706
	707	! Search for index(ni) and size(knon) of domaine to treat
	708	ni(:) = 0
	709	knon = 0
	710	DO i = 1, klon
[996]	711	IF (pctsrf(i,nsrf) > 0.) THEN
[781]	712	knon = knon + 1
	713	ni(knon) = i
	714	ENDIF
	715	ENDDO
	716
	717	! write index, with IOIPSL
	718	IF (debugindex .AND. mpi_size==1) THEN
	719	tabindx(:)=0.
	720	DO i=1,knon
[1403]	721	tabindx(i)=REAL(i)
[781]	722	END DO
	723	debugtab(:,:) = 0.
	724	ndexbg(:) = 0
	725	CALL gath2cpl(tabindx,debugtab,knon,ni)
	726	CALL histwrite(nidbg,cl_surf(nsrf),itap,debugtab,iim*(jjm+1), ndexbg)
	727	ENDIF
	728
	729	!****************************************************************************************
	730	! 5) Compress variables
	731	!
	732	!****************************************************************************************
	733
	734	DO j = 1, knon
	735	i = ni(j)
[888]	736	ypct(j) = pctsrf(i,nsrf)
	737	yts(j) = ts(i,nsrf)
	738	ysnow(j) = snow(i,nsrf)
	739	yqsurf(j) = qsurf(i,nsrf)
	740	yalb(j) = alb(i,nsrf)
	741	yalb1(j) = alb1(i,nsrf)
	742	yalb2(j) = alb2(i,nsrf)
[781]	743	yrain_f(j) = rain_f(i)
	744	ysnow_f(j) = snow_f(i)
	745	yagesno(j) = agesno(i,nsrf)
[888]	746	yfder(j) = fder(i)
	747	ysolsw(j) = solsw(i,nsrf)
	748	ysollw(j) = sollw(i,nsrf)
	749	yrugos(j) = rugos(i,nsrf)
[781]	750	yrugoro(j) = rugoro(i)
[1067]	751	yu1(j) = u(i,1)
	752	yv1(j) = v(i,1)
[781]	753	ypaprs(j,klev+1) = paprs(i,klev+1)
[1067]	754	ywindsp(j) = SQRT(u10m(i,nsrf)2 + v10m(i,nsrf)2 )
[1865]	755	! Martin
	756	yzsig(j) = zsig(i)
	757	ycldt(j) = cldt(i)
	758	yrmu0(j) = rmu0(i)
	759	! Martin
[781]	760	END DO
	761
	762	DO k = 1, klev
	763	DO j = 1, knon
	764	i = ni(j)
	765	ypaprs(j,k) = paprs(i,k)
	766	ypplay(j,k) = pplay(i,k)
[996]	767	ydelp(j,k) = delp(i,k)
	768	ytke(j,k) = tke(i,k,nsrf)
[781]	769	yu(j,k) = u(i,k)
	770	yv(j,k) = v(i,k)
	771	yt(j,k) = t(i,k)
	772	yq(j,k) = q(i,k)
	773	ENDDO
	774	ENDDO
[1761]	775
[781]	776	DO k = 1, nsoilmx
	777	DO j = 1, knon
	778	i = ni(j)
	779	ytsoil(j,k) = ftsoil(i,k,nsrf)
	780	END DO
	781	END DO
	782
	783	! qsol(water height in soil) only for bucket continental model
	784	IF ( nsrf .EQ. is_ter .AND. .NOT. ok_veget ) THEN
	785	DO j = 1, knon
	786	i = ni(j)
	787	yqsol(j) = qsol(i)
	788	END DO
	789	ENDIF
	790
	791	!****************************************************************************************
[1067]	792	! 6a) Calculate coefficients for turbulent diffusion at surface, cdragh et cdragm.
[781]	793	!
	794	!****************************************************************************************
	795
[1067]	796	CALL clcdrag( knon, nsrf, ypaprs, ypplay, &
	797	yu(:,1), yv(:,1), yt(:,1), yq(:,1), &
	798	yts, yqsurf, yrugos, &
	799	ycdragm, ycdragh )
[2126]	800	! --- special Dice: on force cdragm ( a defaut de forcer ustar) MPL 05082013
	801	IF (ok_prescr_ust) then
	802	DO i = 1, knon
	803	print *,'ycdragm avant=',ycdragm(i)
	804	vent= sqrt(yu(i,1)yu(i,1)+yv(i,1)yv(i,1))
	805	! ycdragm(i) = ustust/(1.+(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1)))
	806	! ycdragm(i) = ustust/((1.+sqrt(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1))) &
	807	! sqrt(yu(i,1)yu(i,1)+yv(i,1)*yv(i,1)))
	808	ycdragm(i) = ust*ust/(1.+vent)/vent
	809	print *,'ycdragm ust yu yv apres=',ycdragm(i),ust,yu(i,1),yv(i,1)
	810	ENDDO
	811	ENDIF
[1067]	812
[2126]	813
[1067]	814	!****************************************************************************************
	815	! 6b) Calculate coefficients for turbulent diffusion in the atmosphere, ycoefm et ycoefm.
	816	!
	817	!****************************************************************************************
	818
[781]	819	CALL coef_diff_turb(dtime, nsrf, knon, ni, &
[1067]	820	ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, ycdragm, &
[996]	821	ycoefm, ycoefh, ytke)
[1761]	822
	823	IF (iflag_pbl>=20.AND.iflag_pbl<30) THEN
	824	! In this case, coef_diff_turb is called for the Cd only
	825	DO k = 2, klev
	826	DO j = 1, knon
	827	i = ni(j)
	828	ycoefh(j,k) = zcoefh(i,k,nsrf)
	829	ycoefm(j,k) = zcoefm(i,k,nsrf)
	830	ENDDO
	831	ENDDO
	832	ENDIF
[781]	833
	834	!****************************************************************************************
	835	!
	836	! 8) "La descente" - "The downhill"
	837	!
	838	! climb_hq_down and climb_wind_down calculate the coefficients
	839	! Ccoef_X et Dcoef_X for X=[H, Q, U, V].
	840	! Only the coefficients at surface for H and Q are returned.
	841	!
	842	!****************************************************************************************
	843
	844	! - Calculate the coefficients Ccoef_H, Ccoef_Q, Dcoef_H and Dcoef_Q
	845	CALL climb_hq_down(knon, ycoefh, ypaprs, ypplay, &
	846	ydelp, yt, yq, dtime, &
[1067]	847	AcoefH, AcoefQ, BcoefH, BcoefQ)
[781]	848
	849	! - Calculate the coefficients Ccoef_U, Ccoef_V, Dcoef_U and Dcoef_V
[1067]	850	CALL climb_wind_down(knon, dtime, ycoefm, ypplay, ypaprs, yt, ydelp, yu, yv, &
	851	AcoefU, AcoefV, BcoefU, BcoefV)
[781]	852
	853
	854	!****************************************************************************************
	855	! 9) Small calculations
	856	!
	857	!****************************************************************************************
[888]	858
	859	! - Reference pressure is given the values at surface level
[781]	860	ypsref(:) = ypaprs(:,1)
	861
[1279]	862	! - CO2 field on 2D grid to be sent to ORCHIDEE
	863	! Transform to compressed field
	864	IF (carbon_cycle_cpl) THEN
	865	DO i=1,knon
	866	r_co2_ppm(i) = co2_send(ni(i))
	867	END DO
	868	ELSE
	869	r_co2_ppm(:) = co2_ppm ! Constant field
	870	END IF
[781]	871
	872	!****************************************************************************************
	873	!
[1146]	874	! Calulate t2m and q2m for the case of calculation at land grid points
	875	! t2m and q2m are needed as input to ORCHIDEE
	876	!
	877	!****************************************************************************************
	878	IF (nsrf == is_ter) THEN
	879
	880	DO i = 1, knon
	881	zgeo1(i) = RD * yt(i,1) / (0.5*(ypaprs(i,1)+ypplay(i,1))) &
	882	* (ypaprs(i,1)-ypplay(i,1))
	883	END DO
	884
	885	! Calculate the temperature et relative humidity at 2m and the wind at 10m
	886	CALL stdlevvar(klon, knon, is_ter, zxli, &
	887	yu(:,1), yv(:,1), yt(:,1), yq(:,1), zgeo1, &
	888	yts, yqsurf, yrugos, ypaprs(:,1), ypplay(:,1), &
	889	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
	890
	891	END IF
	892
	893	!****************************************************************************************
	894	!
[781]	895	! 10) Switch selon current surface
	896	! It is necessary to start with the continental surfaces because the ocean
	897	! needs their run-off.
	898	!
	899	!****************************************************************************************
	900	SELECT CASE(nsrf)
	901
	902	CASE(is_ter)
[888]	903	! ylwdown : to be removed, calculation is now done at land surface in surf_land
	904	ylwdown(:)=0.0
	905	DO i=1,knon
	906	ylwdown(i)=lwdown_m(ni(i))
	907	END DO
[781]	908	CALL surf_land(itap, dtime, date0, jour, knon, ni,&
	909	rlon, rlat, &
[888]	910	debut, lafin, ydelp(:,1), r_co2_ppm, ysolsw, ysollw, yalb, &
[1067]	911	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
	912	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	913	AcoefU, AcoefV, BcoefU, BcoefV, &
[781]	914	ypsref, yu1, yv1, yrugoro, pctsrf, &
[1146]	915	ylwdown, yq2m, yt2m, &
[888]	916	ysnow, yqsol, yagesno, ytsoil, &
	917	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
[996]	918	yqsurf, ytsurf_new, y_dflux_t, y_dflux_q, &
[1146]	919	y_flux_u1, y_flux_v1 )
	920
[2126]	921	! Special DICE MPL 05082013
	922	IF (ok_prescr_ust) THEN
	923	! ysnow(:)=0.
	924	! yqsol(:)=0.
	925	! yagesno(:)=50.
	926	! ytsoil(:,:)=300.
	927	! yz0_new(:)=0.001
	928	! yalb1_new(:)=0.22
	929	! yalb2_new(:)=0.22
	930	! yevap(:)=flat/RLVTT
	931	! yfluxlat(:)=-flat
	932	! yfluxsens(:)=-fsens
	933	! yqsurf(:)=0.
	934	! ytsurf_new(:)=tg
	935	! y_dflux_t(:)=0.
	936	! y_dflux_q(:)=0.
	937	y_flux_u1(:)=ycdragm(:)(1.+sqrt(yu(:,1)yu(:,1)+yv(:,1)yv(:,1)))yu(:,1)*ypplay(:,1)/RD/yt(:,1)
	938	y_flux_v1(:)=ycdragm(:)(1.+sqrt(yu(:,1)yu(:,1)+yv(:,1)yv(:,1)))yv(:,1)*ypplay(:,1)/RD/yt(:,1)
	939	ENDIF
	940
[781]	941
	942	CASE(is_lic)
[1865]	943	! Martin
[781]	944	CALL surf_landice(itap, dtime, knon, ni, &
[1865]	945	rlon, rlat, debut, lafin, &
	946	yrmu0, ysollwd, yalb, ypphi(:,1), &
[888]	947	ysolsw, ysollw, yts, ypplay(:,1), &
[1067]	948	ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
	949	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	950	AcoefU, AcoefV, BcoefU, BcoefV, &
[781]	951	ypsref, yu1, yv1, yrugoro, pctsrf, &
[888]	952	ysnow, yqsurf, yqsol, yagesno, &
	953	ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
[1067]	954	ytsurf_new, y_dflux_t, y_dflux_q, &
[1865]	955	yzsig, ycldt, &
	956	ysnowhgt, yqsnow, ytoice, ysissnow, &
	957	yalb3_new, yrunoff, &
[1067]	958	y_flux_u1, y_flux_v1)
[1865]	959	!CALL surf_landice(itap, dtime, knon, ni, &
	960	! ysolsw, ysollw, yts, ypplay(:,1), &
	961	! ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
	962	! AcoefH, AcoefQ, BcoefH, BcoefQ, &
	963	! AcoefU, AcoefV, BcoefU, BcoefV, &
	964	! ypsref, yu1, yv1, yrugoro, pctsrf, &
	965	! ysnow, yqsurf, yqsol, yagesno, &
	966	! ytsoil, yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
	967	! ytsurf_new, y_dflux_t, y_dflux_q, &
	968	! y_flux_u1, y_flux_v1)
	969
[1872]	970	alb3_lic(:)=0.
[1865]	971	DO j = 1, knon
	972	i = ni(j)
	973	alb3_lic(i) = yalb3_new(j)
	974	snowhgt(i) = ysnowhgt(j)
	975	qsnow(i) = yqsnow(j)
	976	to_ice(i) = ytoice(j)
	977	sissnow(i) = ysissnow(j)
	978	runoff(i) = yrunoff(j)
	979	END DO
	980	! Martin
[781]	981
	982	CASE(is_oce)
[888]	983	CALL surf_ocean(rlon, rlat, ysolsw, ysollw, yalb1, &
[996]	984	yrugos, ywindsp, rmu0, yfder, yts, &
[781]	985	itap, dtime, jour, knon, ni, &
[1067]	986	ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
	987	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	988	AcoefU, AcoefV, BcoefU, BcoefV, &
[781]	989	ypsref, yu1, yv1, yrugoro, pctsrf, &
[888]	990	ysnow, yqsurf, yagesno, &
	991	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
[1067]	992	ytsurf_new, y_dflux_t, y_dflux_q, slab_wfbils, &
	993	y_flux_u1, y_flux_v1)
[781]	994
	995	CASE(is_sic)
	996	CALL surf_seaice( &
[888]	997	rlon, rlat, ysolsw, ysollw, yalb1, yfder, &
[781]	998	itap, dtime, jour, knon, ni, &
[1067]	999	lafin, &
	1000	yts, ypplay(:,1), ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&
	1001	AcoefH, AcoefQ, BcoefH, BcoefQ, &
	1002	AcoefU, AcoefV, BcoefU, BcoefV, &
[781]	1003	ypsref, yu1, yv1, yrugoro, pctsrf, &
[888]	1004	ysnow, yqsurf, yqsol, yagesno, ytsoil, &
	1005	yz0_new, yalb1_new, yalb2_new, yevap, yfluxsens, yfluxlat, &
[1067]	1006	ytsurf_new, y_dflux_t, y_dflux_q, &
	1007	y_flux_u1, y_flux_v1)
[781]	1008
	1009
	1010	CASE DEFAULT
	1011	WRITE(lunout,*) 'Surface index = ', nsrf
	1012	abort_message = 'Surface index not valid'
	1013	CALL abort_gcm(modname,abort_message,1)
	1014	END SELECT
	1015
	1016
	1017	!****************************************************************************************
	1018	! 11) - Calcul the increment of surface temperature
	1019	!
	1020	!****************************************************************************************
[1894]	1021
	1022	if (evap0>=0.) then
	1023	yevap(:)=evap0
	1024	yevap(:)=RLVTT*evap0
	1025	endif
	1026
	1027
[781]	1028	y_d_ts(1:knon) = ytsurf_new(1:knon) - yts(1:knon)
	1029
	1030	!****************************************************************************************
	1031	!
	1032	! 12) "La remontee" - "The uphill"
	1033	!
	1034	! The fluxes (y_flux_X) and tendancy (y_d_X) are calculated
	1035	! for X=H, Q, U and V, for all vertical levels.
	1036	!
	1037	!****************************************************************************************
	1038	! H and Q
[1067]	1039	IF (ok_flux_surf) THEN
	1040	PRINT *,'pbl_surface: fsens flat RLVTT=',fsens,flat,RLVTT
[882]	1041	y_flux_t1(:) = fsens
	1042	y_flux_q1(:) = flat/RLVTT
	1043	yfluxlat(:) = flat
[1555]	1044
	1045	Kech_h(:) = ycdragh(:) * (1.0+SQRT(yu(:,1)2+yv(:,1)2)) * &
	1046	ypplay(:,1)/(RD*yt(:,1))
	1047	ytoto(:)=(1./RCPD)(AcoefH(:)+BcoefH(:)y_flux_t1(:)*dtime)
	1048	ytsurf_new(:)=ytoto(:)-y_flux_t1(:)/(Kech_h(:)*RCPD)
	1049	y_d_ts(:) = ytsurf_new(:) - yts(:)
	1050
[1067]	1051	ELSE
[882]	1052	y_flux_t1(:) = yfluxsens(:)
	1053	y_flux_q1(:) = -yevap(:)
[1067]	1054	ENDIF
[781]	1055
	1056	CALL climb_hq_up(knon, dtime, yt, yq, &
	1057	y_flux_q1, y_flux_t1, ypaprs, ypplay, &
	1058	y_flux_q(:,:), y_flux_t(:,:), y_d_q(:,:), y_d_t(:,:))
	1059
[1067]	1060
	1061	CALL climb_wind_up(knon, dtime, yu, yv, y_flux_u1, y_flux_v1, &
[781]	1062	y_flux_u, y_flux_v, y_d_u, y_d_v)
	1063
[1067]	1064
[1761]	1065	y_d_t_diss(:,:)=0.
	1066	IF (iflag_pbl>=20 .and. iflag_pbl<30) THEN
	1067	CALL yamada_c(knon,dtime,ypaprs,ypplay &
	1068	& ,yu,yv,yt,y_d_u,y_d_v,y_d_t,ycdragm,ytke,ycoefm,ycoefh,ycoefq,y_d_t_diss,yustar &
	1069	& ,iflag_pbl,nsrf)
	1070	ENDIF
	1071	! print*,'yamada_c OK'
	1072
[781]	1073	DO j = 1, knon
	1074	y_dflux_t(j) = y_dflux_t(j) * ypct(j)
	1075	y_dflux_q(j) = y_dflux_q(j) * ypct(j)
	1076	ENDDO
	1077
	1078	!****************************************************************************************
	1079	! 13) Transform variables for output format :
	1080	! - Decompress
	1081	! - Multiply with pourcentage of current surface
	1082	! - Cumulate in global variable
	1083	!
	1084	!****************************************************************************************
	1085
	1086	DO k = 1, klev
	1087	DO j = 1, knon
	1088	i = ni(j)
[1761]	1089	y_d_t_diss(j,k) = y_d_t_diss(j,k) * ypct(j)
[996]	1090	y_d_t(j,k) = y_d_t(j,k) * ypct(j)
	1091	y_d_q(j,k) = y_d_q(j,k) * ypct(j)
	1092	y_d_u(j,k) = y_d_u(j,k) * ypct(j)
	1093	y_d_v(j,k) = y_d_v(j,k) * ypct(j)
[781]	1094
	1095	flux_t(i,k,nsrf) = y_flux_t(j,k)
	1096	flux_q(i,k,nsrf) = y_flux_q(j,k)
	1097	flux_u(i,k,nsrf) = y_flux_u(j,k)
	1098	flux_v(i,k,nsrf) = y_flux_v(j,k)
[878]	1099
	1100
[781]	1101	ENDDO
	1102	ENDDO
[1067]	1103
[1761]	1104	! print*,'Dans pbl OK1'
	1105
[781]	1106	evap(:,nsrf) = - flux_q(:,1,nsrf)
	1107
[888]	1108	alb1(:, nsrf) = 0.
	1109	alb2(:, nsrf) = 0.
[781]	1110	snow(:, nsrf) = 0.
	1111	qsurf(:, nsrf) = 0.
	1112	rugos(:, nsrf) = 0.
	1113	fluxlat(:,nsrf) = 0.
	1114	DO j = 1, knon
	1115	i = ni(j)
	1116	d_ts(i,nsrf) = y_d_ts(j)
[888]	1117	alb1(i,nsrf) = yalb1_new(j)
	1118	alb2(i,nsrf) = yalb2_new(j)
[781]	1119	snow(i,nsrf) = ysnow(j)
	1120	qsurf(i,nsrf) = yqsurf(j)
	1121	rugos(i,nsrf) = yz0_new(j)
	1122	fluxlat(i,nsrf) = yfluxlat(j)
	1123	agesno(i,nsrf) = yagesno(j)
[1067]	1124	cdragh(i) = cdragh(i) + ycdragh(j)*ypct(j)
	1125	cdragm(i) = cdragm(i) + ycdragm(j)*ypct(j)
[781]	1126	dflux_t(i) = dflux_t(i) + y_dflux_t(j)
	1127	dflux_q(i) = dflux_q(i) + y_dflux_q(j)
	1128	END DO
	1129
[1761]	1130	! print*,'Dans pbl OK2'
	1131
[1067]	1132	DO k = 2, klev
	1133	DO j = 1, knon
	1134	i = ni(j)
[1761]	1135	tke(i,k,nsrf) = ytke(j,k)
	1136	zcoefh(i,k,nsrf) = ycoefh(j,k)
	1137	zcoefm(i,k,nsrf) = ycoefm(j,k)
	1138	tke(i,k,is_ave) = tke(i,k,is_ave) + ytke(j,k)*ypct(j)
	1139	zcoefh(i,k,is_ave) = zcoefh(i,k,is_ave) + ycoefh(j,k)*ypct(j)
	1140	zcoefm(i,k,is_ave) = zcoefm(i,k,is_ave) + ycoefm(j,k)*ypct(j)
[1067]	1141	END DO
	1142	END DO
	1143
[1761]	1144	! print*,'Dans pbl OK3'
	1145
[781]	1146	IF ( nsrf .EQ. is_ter ) THEN
	1147	DO j = 1, knon
	1148	i = ni(j)
	1149	qsol(i) = yqsol(j)
	1150	END DO
	1151	END IF
	1152
	1153	ftsoil(:,:,nsrf) = 0.
	1154	DO k = 1, nsoilmx
	1155	DO j = 1, knon
	1156	i = ni(j)
	1157	ftsoil(i, k, nsrf) = ytsoil(j,k)
	1158	END DO
	1159	END DO
	1160
	1161
	1162	DO k = 1, klev
	1163	DO j = 1, knon
	1164	i = ni(j)
[1761]	1165	d_t_diss(i,k) = d_t_diss(i,k) + y_d_t_diss(j,k)
[781]	1166	d_t(i,k) = d_t(i,k) + y_d_t(j,k)
	1167	d_q(i,k) = d_q(i,k) + y_d_q(j,k)
	1168	d_u(i,k) = d_u(i,k) + y_d_u(j,k)
	1169	d_v(i,k) = d_v(i,k) + y_d_v(j,k)
	1170	END DO
	1171	END DO
	1172
[1761]	1173	! print*,'Dans pbl OK4'
	1174
[781]	1175	!****************************************************************************************
	1176	! 14) Calculate the temperature et relative humidity at 2m and the wind at 10m
	1177	! Call HBTM
	1178	!
	1179	!****************************************************************************************
	1180	t2m(:,nsrf) = 0.
	1181	q2m(:,nsrf) = 0.
[1670]	1182	ustar(:,nsrf) = 0.
[1816]	1183	wstar(:,nsrf) = 0.
[781]	1184	u10m(:,nsrf) = 0.
	1185	v10m(:,nsrf) = 0.
	1186	pblh(:,nsrf) = 0. ! Hauteur de couche limite
	1187	plcl(:,nsrf) = 0. ! Niveau de condensation de la CLA
	1188	capCL(:,nsrf) = 0. ! CAPE de couche limite
	1189	oliqCL(:,nsrf) = 0. ! eau_liqu integree de couche limite
	1190	cteiCL(:,nsrf) = 0. ! cloud top instab. crit. couche limite
	1191	pblt(:,nsrf) = 0. ! T a la Hauteur de couche limite
	1192	therm(:,nsrf) = 0.
	1193	trmb1(:,nsrf) = 0. ! deep_cape
	1194	trmb2(:,nsrf) = 0. ! inhibition
	1195	trmb3(:,nsrf) = 0. ! Point Omega
	1196
	1197	#undef T2m
	1198	#define T2m
	1199	#ifdef T2m
[996]	1200	! Calculations of diagnostic t,q at 2m and u, v at 10m
[781]	1201
[1761]	1202	! print*,'Dans pbl OK41'
	1203	! print*,'tair1,yt(:,1),y_d_t(:,1)'
	1204	! print*, tair1,yt(:,1),y_d_t(:,1)
[781]	1205	DO j=1, knon
	1206	i = ni(j)
	1207	uzon(j) = yu(j,1) + y_d_u(j,1)
	1208	vmer(j) = yv(j,1) + y_d_v(j,1)
[1761]	1209	tair1(j) = yt(j,1) + y_d_t(j,1) + y_d_t_diss(j,1)
[781]	1210	qair1(j) = yq(j,1) + y_d_q(j,1)
	1211	zgeo1(j) = RD * tair1(j) / (0.5*(ypaprs(j,1)+ypplay(j,1))) &
	1212	* (ypaprs(j,1)-ypplay(j,1))
	1213	tairsol(j) = yts(j) + y_d_ts(j)
	1214	rugo1(j) = yrugos(j)
	1215	IF(nsrf.EQ.is_oce) THEN
	1216	rugo1(j) = rugos(i,nsrf)
	1217	ENDIF
	1218	psfce(j)=ypaprs(j,1)
	1219	patm(j)=ypplay(j,1)
	1220	qairsol(j) = yqsurf(j)
	1221	END DO
	1222
[1761]	1223	! print*,'Dans pbl OK42A'
	1224	! print*,'tair1,yt(:,1),y_d_t(:,1)'
	1225	! print*, tair1,yt(:,1),y_d_t(:,1)
[781]	1226
	1227	! Calculate the temperature et relative humidity at 2m and the wind at 10m
	1228	CALL stdlevvar(klon, knon, nsrf, zxli, &
	1229	uzon, vmer, tair1, qair1, zgeo1, &
	1230	tairsol, qairsol, rugo1, psfce, patm, &
	1231	yt2m, yq2m, yt10m, yq10m, yu10m, yustar)
[1761]	1232	! print*,'Dans pbl OK42B'
[781]	1233
	1234	DO j=1, knon
	1235	i = ni(j)
	1236	t2m(i,nsrf)=yt2m(j)
[996]	1237	q2m(i,nsrf)=yq2m(j)
[781]	1238
[996]	1239	! u10m, v10m : composantes du vent a 10m sans spirale de Ekman
[1670]	1240	ustar(i,nsrf)=yustar(j)
[781]	1241	u10m(i,nsrf)=(yu10m(j) * uzon(j))/SQRT(uzon(j)2+vmer(j)2)
	1242	v10m(i,nsrf)=(yu10m(j) * vmer(j))/SQRT(uzon(j)2+vmer(j)2)
[1670]	1243
[781]	1244	END DO
	1245
[1761]	1246	! print*,'Dans pbl OK43'
[996]	1247	!IM Calcule de l'humidite relative a 2m (rh2m) pour diagnostique
	1248	!IM Ajoute dependance type surface
	1249	IF (thermcep) THEN
	1250	DO j = 1, knon
	1251	i=ni(j)
	1252	zdelta1 = MAX(0.,SIGN(1., rtt-yt2m(j) ))
	1253	zx_qs1 = r2es * FOEEW(yt2m(j),zdelta1)/paprs(i,1)
	1254	zx_qs1 = MIN(0.5,zx_qs1)
	1255	zcor1 = 1./(1.-RETV*zx_qs1)
	1256	zx_qs1 = zx_qs1*zcor1
	1257
	1258	rh2m(i) = rh2m(i) + yq2m(j)/zx_qs1 * pctsrf(i,nsrf)
	1259	qsat2m(i) = qsat2m(i) + zx_qs1 * pctsrf(i,nsrf)
	1260	END DO
	1261	END IF
[781]	1262
[1761]	1263	! print*,'OK pbl 5'
[1816]	1264	CALL hbtm(knon, ypaprs, ypplay, &
	1265	yt2m,yt10m,yq2m,yq10m,yustar,ywstar, &
[781]	1266	y_flux_t,y_flux_q,yu,yv,yt,yq, &
	1267	ypblh,ycapCL,yoliqCL,ycteiCL,ypblT, &
	1268	ytherm,ytrmb1,ytrmb2,ytrmb3,ylcl)
	1269
	1270	DO j=1, knon
	1271	i = ni(j)
	1272	pblh(i,nsrf) = ypblh(j)
[1816]	1273	wstar(i,nsrf) = ywstar(j)
[781]	1274	plcl(i,nsrf) = ylcl(j)
	1275	capCL(i,nsrf) = ycapCL(j)
	1276	oliqCL(i,nsrf) = yoliqCL(j)
	1277	cteiCL(i,nsrf) = ycteiCL(j)
	1278	pblT(i,nsrf) = ypblT(j)
	1279	therm(i,nsrf) = ytherm(j)
	1280	trmb1(i,nsrf) = ytrmb1(j)
	1281	trmb2(i,nsrf) = ytrmb2(j)
	1282	trmb3(i,nsrf) = ytrmb3(j)
	1283	END DO
	1284
[1761]	1285	! print*,'OK pbl 6'
[781]	1286	#else
[996]	1287	! T2m not defined
[781]	1288	! No calculation
[996]	1289	PRINT*,' Warning !!! No T2m calculation. Output is set to zero.'
[781]	1290	#endif
	1291
	1292	!****************************************************************************************
	1293	! 15) End of loop over different surfaces
	1294	!
	1295	!****************************************************************************************
	1296	END DO loop_nbsrf
	1297
	1298	!****************************************************************************************
	1299	! 16) Calculate the mean value over all sub-surfaces for som variables
	1300	!
	1301	!****************************************************************************************
	1302
[1761]	1303	! print*,'OK pbl 7'
[781]	1304	zxfluxt(:,:) = 0.0 ; zxfluxq(:,:) = 0.0
	1305	zxfluxu(:,:) = 0.0 ; zxfluxv(:,:) = 0.0
	1306	DO nsrf = 1, nbsrf
	1307	DO k = 1, klev
	1308	DO i = 1, klon
[996]	1309	zxfluxt(i,k) = zxfluxt(i,k) + flux_t(i,k,nsrf) * pctsrf(i,nsrf)
	1310	zxfluxq(i,k) = zxfluxq(i,k) + flux_q(i,k,nsrf) * pctsrf(i,nsrf)
	1311	zxfluxu(i,k) = zxfluxu(i,k) + flux_u(i,k,nsrf) * pctsrf(i,nsrf)
	1312	zxfluxv(i,k) = zxfluxv(i,k) + flux_v(i,k,nsrf) * pctsrf(i,nsrf)
[781]	1313	END DO
	1314	END DO
	1315	END DO
	1316
[1761]	1317	! print*,'OK pbl 8'
[781]	1318	DO i = 1, klon
	1319	zxsens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
	1320	zxevap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
	1321	fder_print(i) = fder(i) + dflux_t(i) + dflux_q(i)
	1322	ENDDO
	1323
	1324	!
	1325	! Incrementer la temperature du sol
	1326	!
	1327	zxtsol(:) = 0.0 ; zxfluxlat(:) = 0.0
	1328	zt2m(:) = 0.0 ; zq2m(:) = 0.0
[1670]	1329	zustar(:)=0.0 ; zu10m(:) = 0.0 ; zv10m(:) = 0.0
[781]	1330	s_pblh(:) = 0.0 ; s_plcl(:) = 0.0
	1331	s_capCL(:) = 0.0 ; s_oliqCL(:) = 0.0
	1332	s_cteiCL(:) = 0.0; s_pblT(:) = 0.0
	1333	s_therm(:) = 0.0 ; s_trmb1(:) = 0.0
	1334	s_trmb2(:) = 0.0 ; s_trmb3(:) = 0.0
[1816]	1335	wstar(:,is_ave)=0.
[781]	1336
[1761]	1337	! print*,'OK pbl 9'
[781]	1338
	1339	DO nsrf = 1, nbsrf
	1340	DO i = 1, klon
	1341	ts(i,nsrf) = ts(i,nsrf) + d_ts(i,nsrf)
	1342
	1343	wfbils(i,nsrf) = ( solsw(i,nsrf) + sollw(i,nsrf) &
[996]	1344	+ flux_t(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
[781]	1345	wfbilo(i,nsrf) = (evap(i,nsrf) - (rain_f(i) + snow_f(i))) * &
[996]	1346	pctsrf(i,nsrf)
[781]	1347
[996]	1348	zxtsol(i) = zxtsol(i) + ts(i,nsrf) * pctsrf(i,nsrf)
	1349	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf) * pctsrf(i,nsrf)
[781]	1350
[996]	1351	zt2m(i) = zt2m(i) + t2m(i,nsrf) * pctsrf(i,nsrf)
	1352	zq2m(i) = zq2m(i) + q2m(i,nsrf) * pctsrf(i,nsrf)
[1670]	1353	zustar(i) = zustar(i) + ustar(i,nsrf) * pctsrf(i,nsrf)
[1816]	1354	wstar(i,is_ave)=wstar(i,is_ave)+wstar(i,nsrf)*pctsrf(i,nsrf)
[996]	1355	zu10m(i) = zu10m(i) + u10m(i,nsrf) * pctsrf(i,nsrf)
	1356	zv10m(i) = zv10m(i) + v10m(i,nsrf) * pctsrf(i,nsrf)
[781]	1357
[996]	1358	s_pblh(i) = s_pblh(i) + pblh(i,nsrf) * pctsrf(i,nsrf)
	1359	s_plcl(i) = s_plcl(i) + plcl(i,nsrf) * pctsrf(i,nsrf)
	1360	s_capCL(i) = s_capCL(i) + capCL(i,nsrf) * pctsrf(i,nsrf)
	1361	s_oliqCL(i) = s_oliqCL(i) + oliqCL(i,nsrf)* pctsrf(i,nsrf)
	1362	s_cteiCL(i) = s_cteiCL(i) + cteiCL(i,nsrf)* pctsrf(i,nsrf)
	1363	s_pblT(i) = s_pblT(i) + pblT(i,nsrf) * pctsrf(i,nsrf)
	1364	s_therm(i) = s_therm(i) + therm(i,nsrf) * pctsrf(i,nsrf)
	1365	s_trmb1(i) = s_trmb1(i) + trmb1(i,nsrf) * pctsrf(i,nsrf)
	1366	s_trmb2(i) = s_trmb2(i) + trmb2(i,nsrf) * pctsrf(i,nsrf)
	1367	s_trmb3(i) = s_trmb3(i) + trmb3(i,nsrf) * pctsrf(i,nsrf)
[781]	1368	END DO
	1369	END DO
[1761]	1370	! print*,'OK pbl 10'
[781]	1371
	1372	IF (check) THEN
	1373	amn=MIN(ts(1,is_ter),1000.)
	1374	amx=MAX(ts(1,is_ter),-1000.)
	1375	DO i=2, klon
	1376	amn=MIN(ts(i,is_ter),amn)
	1377	amx=MAX(ts(i,is_ter),amx)
	1378	ENDDO
	1379	PRINT*,' debut apres d_ts min max ftsol(ts)',itap,amn,amx
	1380	ENDIF
[1067]	1381
	1382	!jg ?
[996]	1383	!!$!
	1384	!!$! If a sub-surface does not exsist for a grid point, the mean value for all
	1385	!!$! sub-surfaces is distributed.
	1386	!!$!
	1387	!!$ DO nsrf = 1, nbsrf
	1388	!!$ DO i = 1, klon
	1389	!!$ IF ((pctsrf_new(i,nsrf) .LT. epsfra) .OR. (t2m(i,nsrf).EQ.0.)) THEN
	1390	!!$ ts(i,nsrf) = zxtsol(i)
	1391	!!$ t2m(i,nsrf) = zt2m(i)
	1392	!!$ q2m(i,nsrf) = zq2m(i)
	1393	!!$ u10m(i,nsrf) = zu10m(i)
	1394	!!$ v10m(i,nsrf) = zv10m(i)
	1395	!!$
	1396	!!$! Les variables qui suivent sont plus utilise, donc peut-etre pas la peine a les mettre ajour
	1397	!!$ pblh(i,nsrf) = s_pblh(i)
	1398	!!$ plcl(i,nsrf) = s_plcl(i)
	1399	!!$ capCL(i,nsrf) = s_capCL(i)
	1400	!!$ oliqCL(i,nsrf) = s_oliqCL(i)
	1401	!!$ cteiCL(i,nsrf) = s_cteiCL(i)
	1402	!!$ pblT(i,nsrf) = s_pblT(i)
	1403	!!$ therm(i,nsrf) = s_therm(i)
	1404	!!$ trmb1(i,nsrf) = s_trmb1(i)
	1405	!!$ trmb2(i,nsrf) = s_trmb2(i)
	1406	!!$ trmb3(i,nsrf) = s_trmb3(i)
	1407	!!$ ENDIF
	1408	!!$ ENDDO
	1409	!!$ ENDDO
[781]	1410
	1411
	1412	DO i = 1, klon
	1413	fder(i) = - 4.0RSIGMAzxtsol(i)**3
	1414	ENDDO
	1415
	1416	zxqsurf(:) = 0.0
	1417	zxsnow(:) = 0.0
	1418	DO nsrf = 1, nbsrf
	1419	DO i = 1, klon
[996]	1420	zxqsurf(i) = zxqsurf(i) + qsurf(i,nsrf) * pctsrf(i,nsrf)
	1421	zxsnow(i) = zxsnow(i) + snow(i,nsrf) * pctsrf(i,nsrf)
[781]	1422	END DO
	1423	END DO
	1424
[1067]	1425	! Premier niveau de vent sortie dans physiq.F
	1426	zu1(:) = u(:,1)
	1427	zv1(:) = v(:,1)
[781]	1428
	1429	! Some of the module declared variables are returned for printing in physiq.F
	1430	qsol_d(:) = qsol(:)
	1431	evap_d(:,:) = evap(:,:)
	1432	rugos_d(:,:) = rugos(:,:)
	1433	agesno_d(:,:) = agesno(:,:)
	1434
	1435
	1436	END SUBROUTINE pbl_surface
	1437	!
	1438	!****************************************************************************************
	1439	!
	1440	SUBROUTINE pbl_surface_final(qsol_rst, fder_rst, snow_rst, qsurf_rst, &
	1441	evap_rst, rugos_rst, agesno_rst, ftsoil_rst)
	1442
[1785]	1443	USE indice_sol_mod
	1444
[781]	1445	INCLUDE "dimsoil.h"
	1446
	1447	! Ouput variables
	1448	!****************************************************************************************
	1449	REAL, DIMENSION(klon), INTENT(OUT) :: qsol_rst
	1450	REAL, DIMENSION(klon), INTENT(OUT) :: fder_rst
	1451	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: snow_rst
	1452	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: qsurf_rst
	1453	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: evap_rst
	1454	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: rugos_rst
	1455	REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: agesno_rst
	1456	REAL, DIMENSION(klon, nsoilmx, nbsrf), INTENT(OUT) :: ftsoil_rst
	1457
	1458
	1459	!****************************************************************************************
	1460	! Return module variables for writing to restart file
	1461	!
	1462	!****************************************************************************************
	1463	qsol_rst(:) = qsol(:)
	1464	fder_rst(:) = fder(:)
	1465	snow_rst(:,:) = snow(:,:)
	1466	qsurf_rst(:,:) = qsurf(:,:)
	1467	evap_rst(:,:) = evap(:,:)
	1468	rugos_rst(:,:) = rugos(:,:)
	1469	agesno_rst(:,:) = agesno(:,:)
	1470	ftsoil_rst(:,:,:) = ftsoil(:,:,:)
	1471
	1472	!****************************************************************************************
	1473	! Deallocate module variables
	1474	!
	1475	!****************************************************************************************
[1413]	1476	! DEALLOCATE(qsol, fder, snow, qsurf, evap, rugos, agesno, ftsoil)
	1477	IF (ALLOCATED(qsol)) DEALLOCATE(qsol)
	1478	IF (ALLOCATED(fder)) DEALLOCATE(fder)
	1479	IF (ALLOCATED(snow)) DEALLOCATE(snow)
	1480	IF (ALLOCATED(qsurf)) DEALLOCATE(qsurf)
	1481	IF (ALLOCATED(evap)) DEALLOCATE(evap)
	1482	IF (ALLOCATED(rugos)) DEALLOCATE(rugos)
	1483	IF (ALLOCATED(agesno)) DEALLOCATE(agesno)
	1484	IF (ALLOCATED(ftsoil)) DEALLOCATE(ftsoil)
[781]	1485
	1486	END SUBROUTINE pbl_surface_final
	1487	!
	1488	!****************************************************************************************
[996]	1489	!
[1670]	1490	SUBROUTINE pbl_surface_newfrac(itime, pctsrf_new, pctsrf_old, tsurf, alb1, alb2, ustar, u10m, v10m, tke)
[996]	1491
	1492	! Give default values where new fraction has appread
	1493
[1785]	1494	USE indice_sol_mod
	1495
[996]	1496	INCLUDE "dimsoil.h"
	1497	INCLUDE "clesphys.h"
[1236]	1498	INCLUDE "compbl.h"
[996]	1499
	1500	! Input variables
	1501	!****************************************************************************************
	1502	INTEGER, INTENT(IN) :: itime
	1503	REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf_new, pctsrf_old
	1504
	1505	! InOutput variables
	1506	!****************************************************************************************
	1507	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: tsurf
	1508	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: alb1, alb2
[1670]	1509	REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: ustar,u10m, v10m
[996]	1510	REAL, DIMENSION(klon,klev+1,nbsrf), INTENT(INOUT) :: tke
	1511
	1512	! Local variables
	1513	!****************************************************************************************
	1514	INTEGER :: nsrf, nsrf_comp1, nsrf_comp2, nsrf_comp3, i
	1515	CHARACTER(len=80) :: abort_message
	1516	CHARACTER(len=20) :: modname = 'pbl_surface_newfrac'
	1517	INTEGER, DIMENSION(nbsrf) :: nfois=0, mfois=0, pfois=0
	1518	!
	1519	! All at once !!
	1520	!****************************************************************************************
	1521
	1522	DO nsrf = 1, nbsrf
	1523	! First decide complement sub-surfaces
	1524	SELECT CASE (nsrf)
	1525	CASE(is_oce)
	1526	nsrf_comp1=is_sic
	1527	nsrf_comp2=is_ter
	1528	nsrf_comp3=is_lic
	1529	CASE(is_sic)
	1530	nsrf_comp1=is_oce
	1531	nsrf_comp2=is_ter
	1532	nsrf_comp3=is_lic
	1533	CASE(is_ter)
	1534	nsrf_comp1=is_lic
	1535	nsrf_comp2=is_oce
	1536	nsrf_comp3=is_sic
	1537	CASE(is_lic)
	1538	nsrf_comp1=is_ter
	1539	nsrf_comp2=is_oce
	1540	nsrf_comp3=is_sic
	1541	END SELECT
	1542
	1543	! Initialize all new fractions
	1544	DO i=1, klon
	1545	IF (pctsrf_new(i,nsrf) > 0. .AND. pctsrf_old(i,nsrf) == 0.) THEN
[1067]	1546
[996]	1547	IF (pctsrf_old(i,nsrf_comp1) > 0.) THEN
	1548	! Use the complement sub-surface, keeping the continents unchanged
	1549	qsurf(i,nsrf) = qsurf(i,nsrf_comp1)
	1550	evap(i,nsrf) = evap(i,nsrf_comp1)
	1551	rugos(i,nsrf) = rugos(i,nsrf_comp1)
	1552	tsurf(i,nsrf) = tsurf(i,nsrf_comp1)
	1553	alb1(i,nsrf) = alb1(i,nsrf_comp1)
	1554	alb2(i,nsrf) = alb2(i,nsrf_comp1)
[1670]	1555	ustar(i,nsrf) = ustar(i,nsrf_comp1)
[996]	1556	u10m(i,nsrf) = u10m(i,nsrf_comp1)
	1557	v10m(i,nsrf) = v10m(i,nsrf_comp1)
[1236]	1558	if (iflag_pbl > 1) then
	1559	tke(i,:,nsrf) = tke(i,:,nsrf_comp1)
	1560	endif
[996]	1561	mfois(nsrf) = mfois(nsrf) + 1
	1562	ELSE
	1563	! The continents have changed. The new fraction receives the mean sum of the existent fractions
	1564	qsurf(i,nsrf) = qsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + qsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
	1565	evap(i,nsrf) = evap(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + evap(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
	1566	rugos(i,nsrf) = rugos(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + rugos(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
	1567	tsurf(i,nsrf) = tsurf(i,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tsurf(i,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
	1568	alb1(i,nsrf) = alb1(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb1(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
	1569	alb2(i,nsrf) = alb2(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + alb2(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
[1670]	1570	ustar(i,nsrf) = ustar(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + ustar(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
[996]	1571	u10m(i,nsrf) = u10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + u10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
	1572	v10m(i,nsrf) = v10m(i,nsrf_comp2) pctsrf_old(i,nsrf_comp2) + v10m(i,nsrf_comp3) pctsrf_old(i,nsrf_comp3)
[1236]	1573	if (iflag_pbl > 1) then
	1574	tke(i,:,nsrf) = tke(i,:,nsrf_comp2)pctsrf_old(i,nsrf_comp2) + tke(i,:,nsrf_comp3)pctsrf_old(i,nsrf_comp3)
	1575	endif
[996]	1576
	1577	! Security abort. This option has never been tested. To test, comment the following line.
	1578	! abort_message='The fraction of the continents have changed!'
	1579	! CALL abort_gcm(modname,abort_message,1)
	1580	nfois(nsrf) = nfois(nsrf) + 1
	1581	END IF
	1582	snow(i,nsrf) = 0.
	1583	agesno(i,nsrf) = 0.
	1584	ftsoil(i,:,nsrf) = tsurf(i,nsrf)
	1585	ELSE
	1586	pfois(nsrf) = pfois(nsrf)+ 1
	1587	END IF
	1588	END DO
	1589
	1590	END DO
	1591
	1592	END SUBROUTINE pbl_surface_newfrac
	1593
[781]	1594	!
[996]	1595	!****************************************************************************************
	1596	!
[781]	1597
	1598	END MODULE pbl_surface_mod

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

Download in other formats:

Original Format