Context Navigation

pbl_surface_mod.F90 @ 1784

Last change on this file since 1784 was 1413, checked in by idelkadi, 14 years ago
Dans fonte_neige_mod.F90 et pbl_surface_mod.F90 : deallocate des variables sous conditions qu'ls soient allous Dans physiq.F : Nettoyage dans physiq.F des anciennes versions des options iflag_cldcon=5 ou 6.
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 58.5 KB

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

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

Download in other formats:

Original Format