Context Navigation

pbl_surface_mod.F90 @ 1085

Last change on this file since 1085 was 1069, checked in by Laurent Fairhead, 16 years ago

Modification pour utiliser albedo des deux intervalles visible et proche
infra rouge pour la repartition sur sous-surface des flux solaire.

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

File size: 56.8 KB

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

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

Download in other formats:

Original Format