Context Navigation

pbl_surface_mod.F90 @ 1006

Last change on this file since 1006 was 996, checked in by lsce, 16 years ago

Modifications liées au calcul des nouveau sous-fractions
Nettoyage de ocean slab : il reste uniquement la version avec glace de mer forcé
Nouveaux variables pour distiguer la version et type d'ocean : type_ocean=force/slab/couple, version_ocean=opa8/nemo pour couplé ou version_ocean=sicOBS pour slab

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

File size: 57.3 KB

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

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

Download in other formats:

Original Format