Context Navigation

pbl_surface_mod.F90 @ 1203

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

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

Download in other formats:

Original Format