Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

pbl_surface_mod.F90

Last change on this file was 1534, checked in by musat, 13 years ago

Ajouts CFMIP2/CMIP5

6eme fichier de sortie "stations" histstn.nc qui necessite 2 fichiers: PARAM/npCFMIP_param.data contenant le nombre de points (120 pour simulations AMIP, 73 pour aqua) PARAM/pointlocations.txt contenat le numero, les coordonnees (lon,lat) et le nom de chaque station
flag LOGICAL dans tous les appels histwrite_phy pour pouvoir sortir le fichier histstn.nc

NB: 1) les flags de type phys_ que l'on met dans le physiq.def_L* pour ajouter plus de sorties

necessitent dorenavant 6 valeurs, la 6eme correspondant au fichier histstn.nc

2) par defaut le fichier histstn.nc ne sort pas; pour le sortir ajouter les lignes suivantes

dans physiq.def_L*

### Type de fichier : global (n) ou stations (y)
phys_out_filestations = n n n n n y

introduction de 2 jeux de flags pour les taux des GES; taux actuels avec suffixes _act, taux futurs avec "_per" avec 2 appels au rayonnement si taux "_per" different des taux "_act" (utiles pour diags. CFMIP 4CO2)
flags "betaCRF" pour calculs CRF pour experiences sensibilite proprietes optiques eau liquide nuageuse avec initialisations par defaut; sinon besoin de fichier beta_crf.data

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

File size: 58.4 KB

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

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

Download in other formats:

Original Format