Context Navigation

pbl_surface_mod.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 58.6 KB

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

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

Download in other formats:

Original Format