Context Navigation

interface_surf.F90 @ 522

Last change on this file since 522 was 522, checked in by lmdzadmin, 20 years ago
Modifications pour la fermeture en eau (fonte des glaciers) JLD, OM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 94.7 KB

Rev	Line
[223]	1	! $Header$
[90]	2
	3	MODULE interface_surf
	4
	5	! Ce module regroupe toutes les routines gerant l'interface entre le modele
	6	! atmospherique et les modeles de surface (sols continentaux, oceans, glaces)
	7	! Les routines sont les suivantes:
	8	!
	9	! interfsurf_*: routines d'aiguillage vers les interfaces avec les
	10	! differents modeles de surface
	11	! interfsol\
	12	! > routines d'interface proprement dite
	13	! interfoce/
	14	!
	15	! interfstart: routine d'initialisation et de lecture de l'etat initial
	16	! "interface"
	17	! interffin : routine d'ecriture de l'etat de redemmarage de l'interface
	18	!
	19	!
	20	! L. Fairhead, LMD, 02/2000
	21
	22	USE ioipsl
	23
	24	IMPLICIT none
	25
	26	PRIVATE
[224]	27	PUBLIC :: interfsurf,interfsurf_hq, gath2cpl
[90]	28
	29	INTERFACE interfsurf
	30	module procedure interfsurf_hq, interfsurf_vent
	31	END INTERFACE
	32
	33	INTERFACE interfoce
	34	module procedure interfoce_cpl, interfoce_slab, interfoce_lim
	35	END INTERFACE
	36
[109]	37	#include "YOMCST.inc"
[112]	38	#include "indicesol.inc"
[90]	39
[109]	40
[112]	41	! run_off ruissellement total
[522]	42	REAL, ALLOCATABLE, DIMENSION(:),SAVE :: run_off, run_off_lic
[177]	43	real, allocatable, dimension(:),save :: coastalflow, riverflow
[290]	44	!!$PB
[522]	45	REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: tmp_rriv, tmp_rcoa,tmp_rlic
[313]	46	!! pour simuler la fonte des glaciers antarctiques
	47	REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: coeff_iceberg
	48	real, save :: surf_maille
	49	real, save :: cte_flux_iceberg = 6.3e7
	50	integer, save :: num_antarctic = 1
[522]	51	REAL, save :: tau_calv
[290]	52	!!$
[90]	53	CONTAINS
	54	!
	55	!############################################################################
	56	!
[205]	57	SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, &
[177]	58	& klon, iim, jjm, nisurf, knon, knindex, pctsrf, &
	59	& rlon, rlat, cufi, cvfi,&
[441]	60	& debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, qsol,&
[177]	61	& zlev, u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
[90]	62	& tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
[177]	63	& precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &
[171]	64	& fder, taux, tauy, rugos, rugoro, &
[441]	65	& albedo, snow, qsurf, &
[105]	66	& tsurf, p1lay, ps, radsol, &
[112]	67	& ocean, npas, nexca, zmasq, &
[90]	68	& evap, fluxsens, fluxlat, dflux_l, dflux_s, &
[281]	69	& tsol_rad, tsurf_new, alb_new, alblw, emis_new, &
[457]	70	!IM cf. JLD & z0_new, pctsrf_new, agesno)
[522]	71	& z0_new, pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0)
[90]	72
	73
	74	! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general
	75	! (sols continentaux, oceans, glaces) pour les fluxs de chaleur et d'humidite.
	76	! En pratique l'interface se fait entre la couche limite du modele
	77	! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
	78	!
	79	!
	80	! L.Fairhead 02/2000
	81	!
	82	! input:
	83	! itime numero du pas de temps
	84	! klon nombre total de points de grille
[98]	85	! iim, jjm nbres de pts de grille
[90]	86	! dtime pas de temps de la physique (en s)
[205]	87	! date0 jour initial
[109]	88	! jour jour dans l'annee en cours,
	89	! rmu0 cosinus de l'angle solaire zenithal
[101]	90	! nexca pas de temps couplage
[90]	91	! nisurf index de la surface a traiter (1 = sol continental)
	92	! knon nombre de points de la surface a traiter
	93	! knindex index des points de la surface a traiter
[98]	94	! pctsrf tableau des pourcentages de surface de chaque maille
[90]	95	! rlon longitudes
	96	! rlat latitudes
[177]	97	! cufi,cvfi resolution des mailles en x et y (m)
[90]	98	! debut logical: 1er appel a la physique
	99	! lafin logical: dernier appel a la physique
	100	! ok_veget logical: appel ou non au schema de surface continental
	101	! (si false calcul simplifie des fluxs sur les continents)
	102	! zlev hauteur de la premiere couche
	103	! u1_lay vitesse u 1ere couche
	104	! v1_lay vitesse v 1ere couche
	105	! temp_air temperature de l'air 1ere couche
	106	! spechum humidite specifique 1ere couche
[177]	107	! epot_air temp potentielle de l'air
[90]	108	! ccanopy concentration CO2 canopee
	109	! tq_cdrag cdrag
	110	! petAcoef coeff. A de la resolution de la CL pour t
	111	! peqAcoef coeff. A de la resolution de la CL pour q
	112	! petBcoef coeff. B de la resolution de la CL pour t
	113	! peqBcoef coeff. B de la resolution de la CL pour q
	114	! precip_rain precipitation liquide
	115	! precip_snow precipitation solide
[177]	116	! sollw flux IR net a la surface
	117	! sollwdown flux IR descendant a la surface
[90]	118	! swnet flux solaire net
	119	! swdown flux solaire entrant a la surface
[98]	120	! albedo albedo de la surface
[90]	121	! tsurf temperature de surface
	122	! p1lay pression 1er niveau (milieu de couche)
	123	! ps pression au sol
	124	! radsol rayonnement net aus sol (LW + SW)
	125	! ocean type d'ocean utilise (force, slab, couple)
[101]	126	! fder derivee des flux (pour le couplage)
	127	! taux, tauy tension de vents
[157]	128	! rugos rugosite
[105]	129	! zmasq masque terre/ocean
[171]	130	! rugoro rugosite orographique
[522]	131	! run_off_lic_0 runoff glacier du pas de temps precedent
[90]	132	!
	133	! output:
	134	! evap evaporation totale
	135	! fluxsens flux de chaleur sensible
	136	! fluxlat flux de chaleur latente
	137	! tsol_rad
	138	! tsurf_new temperature au sol
	139	! alb_new albedo
	140	! emis_new emissivite
	141	! z0_new surface roughness
[98]	142	! pctsrf_new nouvelle repartition des surfaces
[90]	143
	144
	145	! Parametres d'entree
	146	integer, intent(IN) :: itime
[98]	147	integer, intent(IN) :: iim, jjm
[90]	148	integer, intent(IN) :: klon
	149	real, intent(IN) :: dtime
[205]	150	real, intent(IN) :: date0
[90]	151	integer, intent(IN) :: jour
[109]	152	real, intent(IN) :: rmu0(klon)
[90]	153	integer, intent(IN) :: nisurf
	154	integer, intent(IN) :: knon
[147]	155	integer, dimension(klon), intent(in) :: knindex
[98]	156	real, dimension(klon,nbsrf), intent(IN) :: pctsrf
[90]	157	logical, intent(IN) :: debut, lafin, ok_veget
	158	real, dimension(klon), intent(IN) :: rlon, rlat
[177]	159	real, dimension(klon), intent(IN) :: cufi, cvfi
	160	real, dimension(klon), intent(INOUT) :: tq_cdrag
[147]	161	real, dimension(klon), intent(IN) :: zlev
	162	real, dimension(klon), intent(IN) :: u1_lay, v1_lay
	163	real, dimension(klon), intent(IN) :: temp_air, spechum
[177]	164	real, dimension(klon), intent(IN) :: epot_air, ccanopy
	165	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
[147]	166	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
	167	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
[177]	168	real, dimension(klon), intent(IN) :: sollw, sollwdown, swnet, swdown
	169	real, dimension(klon), intent(IN) :: ps, albedo
[457]	170	!IM cf LF
	171	! real, dimension(klon), intent(INOUT) :: tsurf
	172	! real, dimension(klon), intent(IN) :: p1lay
[147]	173	real, dimension(klon), intent(IN) :: tsurf, p1lay
[235]	174	REAL, DIMENSION(klon), INTENT(INOUT) :: radsol,fder
[98]	175	real, dimension(klon), intent(IN) :: zmasq
[235]	176	real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro
[90]	177	character (len = 6) :: ocean
[112]	178	integer :: npas, nexca ! nombre et pas de temps couplage
[441]	179	real, dimension(klon), intent(INOUT) :: evap, snow, qsurf
[177]	180	!! PB ajout pour soil
	181	logical :: soil_model
	182	integer :: nsoilmx
	183	REAL, DIMENSION(klon, nsoilmx) :: tsoil
[441]	184	REAL, dimension(klon), intent(INOUT) :: qsol
[177]	185	REAL, dimension(klon) :: soilcap
	186	REAL, dimension(klon) :: soilflux
[90]	187	! Parametres de sortie
[147]	188	real, dimension(klon), intent(OUT):: fluxsens, fluxlat
	189	real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new
[281]	190	real, dimension(klon), intent(OUT):: alblw
[147]	191	real, dimension(klon), intent(OUT):: emis_new, z0_new
	192	real, dimension(klon), intent(OUT):: dflux_l, dflux_s
[90]	193	real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
[112]	194	real, dimension(klon), intent(INOUT):: agesno
[522]	195	real, dimension(klon), intent(INOUT):: run_off_lic_0
[90]	196
[457]	197	! Flux thermique utiliser pour fondre la neige
	198	!jld a rajouter real, dimension(klon), intent(INOUT):: ffonte
	199	real, dimension(klon), intent(INOUT):: ffonte
	200	! Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
	201	! hauteur de neige, en kg/m2/s
	202	!jld a rajouter real, dimension(klon), intent(INOUT):: fqcalving
	203	real, dimension(klon), intent(INOUT):: fqcalving
	204
[90]	205	! Local
[179]	206	character (len = 20),save :: modname = 'interfsurf_hq'
[90]	207	character (len = 80) :: abort_message
	208	logical, save :: first_call = .true.
[179]	209	integer, save :: error
[236]	210	integer :: ii, index
[458]	211	logical,save :: check = .false.
[147]	212	real, dimension(klon):: cal, beta, dif_grnd, capsol
[177]	213	!!$PB real, parameter :: calice=1.0/(5.1444e+060.15), tau_gl=86400.5.
	214	real, parameter :: calice=1.0/(5.1444e+060.15), tau_gl=86400.5.
[98]	215	real, parameter :: calsno=1./(2.3867e+06*.15)
[147]	216	real, dimension(klon):: alb_ice
	217	real, dimension(klon):: tsurf_temp
[441]	218	real, dimension(klon):: qsurf_new
[258]	219	!! real, allocatable, dimension(:), save :: alb_neig_grid
[157]	220	real, dimension(klon):: alb_neig, alb_eau
[147]	221	real, DIMENSION(klon):: zfra
[177]	222	logical :: cumul = .false.
[458]	223	INTEGER,dimension(1) :: iloc
	224	INTEGER :: isize
	225	real, dimension(klon):: fder_prev
[522]	226	REAL, dimension(klon) :: bidule
[90]	227
	228	if (check) write(,) 'Entree ', modname
	229	!
	230	! On doit commencer par appeler les schemas de surfaces continentales
	231	! car l'ocean a besoin du ruissellement qui est y calcule
	232	!
	233	if (first_call) then
[522]	234	call conf_interface(tau_calv)
[90]	235	if (nisurf /= is_ter .and. klon > 1) then
	236	write(,)' * Warning *'
	237	write(,)' nisurf = ',nisurf,' /= is_ter = ',is_ter
	238	write(,)'or on doit commencer par les surfaces continentales'
	239	abort_message='voir ci-dessus'
	240	call abort_gcm(modname,abort_message,1)
	241	endif
	242	if (ocean /= 'slab ' .and. ocean /= 'force ' .and. ocean /= 'couple') then
	243	write(,)' * Warning *'
	244	write(,)'Option couplage pour l''ocean = ', ocean
	245	abort_message='option pour l''ocean non valable'
	246	call abort_gcm(modname,abort_message,1)
	247	endif
[105]	248	if ( is_oce > is_sic ) then
	249	write(,)' * Warning *'
	250	write(,)' Pour des raisons de sequencement dans le code'
	251	write(,)' l''ocean doit etre traite avant la banquise'
	252	write(,)' or is_oce = ',is_oce, '> is_sic = ',is_sic
	253	abort_message='voir ci-dessus'
	254	call abort_gcm(modname,abort_message,1)
[112]	255	endif
[262]	256	! allocate(alb_neig_grid(klon), stat = error)
	257	! if (error /= 0) then
	258	! abort_message='Pb allocation alb_neig_grid'
	259	! call abort_gcm(modname,abort_message,1)
	260	! endif
[90]	261	endif
	262	first_call = .false.
[98]	263
[157]	264	! Initialisations diverses
	265	!
[235]	266	!!$ cal=0.; beta=1.; dif_grnd=0.; capsol=0.
	267	!!$ alb_new = 0.; z0_new = 0.; alb_neig = 0.0
	268	!!$! PB
	269	!!$ tsurf_new = 0.
[157]	270
[467]	271	!IM cf JLD
	272	ffonte(1:knon)=0.
	273	fqcalving(1:knon)=0.
	274
[235]	275	cal = 999999. ; beta = 999999. ; dif_grnd = 999999. ; capsol = 999999.
	276	alb_new = 999999. ; z0_new = 999999. ; alb_neig = 999999.
	277	tsurf_new = 999999.
[281]	278	alblw = 999999.
[90]	279	! Aiguillage vers les differents schemas de surface
	280
	281	if (nisurf == is_ter) then
	282	!
	283	! Surface "terre" appel a l'interface avec les sols continentaux
	284	!
	285	! allocation du run-off
[177]	286	if (.not. allocated(coastalflow)) then
	287	allocate(coastalflow(knon), stat = error)
	288	if (error /= 0) then
	289	abort_message='Pb allocation coastalflow'
	290	call abort_gcm(modname,abort_message,1)
	291	endif
	292	allocate(riverflow(knon), stat = error)
	293	if (error /= 0) then
	294	abort_message='Pb allocation riverflow'
	295	call abort_gcm(modname,abort_message,1)
	296	endif
[90]	297	allocate(run_off(knon), stat = error)
	298	if (error /= 0) then
[275]	299	abort_message='Pb allocation run_off'
[90]	300	call abort_gcm(modname,abort_message,1)
	301	endif
[290]	302	!!$PB
[398]	303	ALLOCATE (tmp_rriv(iim,jjm+1), stat=error)
	304	if (error /= 0) then
	305	abort_message='Pb allocation tmp_rriv'
	306	call abort_gcm(modname,abort_message,1)
	307	endif
	308	ALLOCATE (tmp_rcoa(iim,jjm+1), stat=error)
	309	if (error /= 0) then
	310	abort_message='Pb allocation tmp_rcoa'
	311	call abort_gcm(modname,abort_message,1)
	312	endif
[522]	313	ALLOCATE (tmp_rlic(iim,jjm+1), stat=error)
	314	if (error /= 0) then
	315	abort_message='Pb allocation tmp_rlic'
	316	call abort_gcm(modname,abort_message,1)
	317	endif
	318
[290]	319	!!$
[177]	320	else if (size(coastalflow) /= knon) then
[90]	321	write(,)'Bizarre, le nombre de points continentaux'
[177]	322	write(,)'a change entre deux appels. J''arrete ...'
	323	abort_message='voir ci-dessus'
	324	call abort_gcm(modname,abort_message,1)
[90]	325	endif
[275]	326	coastalflow = 0.
	327	riverflow = 0.
[90]	328	!
[109]	329	! Calcul age de la neige
	330	!
[177]	331	!!$ PB ATTENTION changement ordre des appels
[258]	332	!!$ CALL albsno(klon,agesno,alb_neig_grid)
[235]	333
[98]	334	if (.not. ok_veget) then
	335	!
	336	! calcul albedo: lecture albedo fichier CL puis ajout albedo neige
	337	!
[109]	338	call interfsur_lim(itime, dtime, jour, &
	339	& klon, nisurf, knon, knindex, debut, &
[121]	340	& alb_new, z0_new)
[236]	341	!
[441]	342	! calcul snow et qsurf, hydrol adapté
[177]	343	!
[236]	344	CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
	345
[177]	346	IF (soil_model) THEN
[236]	347	CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)
	348	cal(1:knon) = RCPD / soilcap(1:knon)
	349	radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
[177]	350	ELSE
	351	cal = RCPD * capsol
	352	!!$ cal = capsol
	353	ENDIF
	354	CALL calcul_fluxs( klon, knon, nisurf, dtime, &
	355	& tsurf, p1lay, cal, beta, tq_cdrag, ps, &
[441]	356	& precip_rain, precip_snow, snow, qsurf, &
[177]	357	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
	358	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
	359	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
	360
[181]	361	CALL fonte_neige( klon, knon, nisurf, dtime, &
	362	& tsurf, p1lay, cal, beta, tq_cdrag, ps, &
	363	& precip_rain, precip_snow, snow, qsol, &
	364	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
	365	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
[457]	366	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
	367	!IM cf JLD
[522]	368	& fqcalving,ffonte, run_off_lic_0)
[181]	369
[236]	370
[258]	371	call albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
	372	where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
[295]	373	zfra(1:knon) = max(0.0,min(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
[281]	374	alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
	375	& alb_new(1 : knon)*(1.0-zfra(1:knon))
[258]	376	z0_new = sqrt(z0_new2+rugoro2)
[281]	377	alblw(1 : knon) = alb_new(1 : knon)
[258]	378
[98]	379	else
[258]	380	!! CALL albsno(klon,agesno,alb_neig_grid)
[98]	381	!
	382	! appel a sechiba
	383	!
[233]	384	call interfsol(itime, klon, dtime, date0, nisurf, knon, &
[177]	385	& knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
[90]	386	& debut, lafin, ok_veget, &
[177]	387	& zlev, u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
[90]	388	& tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
[177]	389	& precip_rain, precip_snow, sollwdown, swnet, swdown, &
[360]	390	& tsurf, p1lay/100., ps/100., radsol, &
[90]	391	& evap, fluxsens, fluxlat, &
[281]	392	& tsol_rad, tsurf_new, alb_new, alblw, &
[441]	393	& emis_new, z0_new, dflux_l, dflux_s, qsurf_new)
[233]	394
	395	!
	396	! ajout de la contribution du relief
	397	!
	398	z0_new = SQRT(z0_new2+rugoro2)
[365]	399	!
	400	! mise a jour de l'humidite saturante calculee par ORCHIDEE
[441]	401	qsurf(1:knon) = qsurf_new(1:knon)
[233]	402
[98]	403	endif
[90]	404	!
[139]	405	! Remplissage des pourcentages de surface
	406	!
	407	pctsrf_new(:,nisurf) = pctsrf(:,nisurf)
	408
[90]	409	else if (nisurf == is_oce) then
	410
	411	if (check) write(,)'ocean, nisurf = ',nisurf
[98]	412
	413
[90]	414	!
	415	! Surface "ocean" appel a l'interface avec l'ocean
	416	!
[101]	417	if (ocean == 'couple') then
	418	if (nexca == 0) then
	419	abort_message='nexca = 0 dans interfoce_cpl'
	420	call abort_gcm(modname,abort_message,1)
	421	endif
	422
[177]	423	cumul = .false.
	424
[458]	425	iloc = maxloc(fder(1:klon))
	426	if (check) then
	427	if (fder(iloc(1))> 0.) then
	428	WRITE(,)'** Debug fder **'
	429	WRITE(,)'max fder(',iloc(1),') = ',fder(iloc(1))
	430	endif
	431	endif
	432	!!$
	433	!!$ where(fder.gt.0.)
	434	!!$ fder = 0.
	435	!!$ endwhere
	436
[177]	437	call interfoce(itime, dtime, cumul, &
[101]	438	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
[112]	439	& ocean, npas, nexca, debut, lafin, &
[177]	440	& swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
	441	& fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
[274]	442	& tsurf_new, alb_new, pctsrf_new)
[101]	443
[90]	444	! else if (ocean == 'slab ') then
	445	! call interfoce(nisurf)
[109]	446	else ! lecture conditions limites
	447	call interfoce(itime, dtime, jour, &
	448	& klon, nisurf, knon, knindex, &
	449	& debut, &
	450	& tsurf_new, pctsrf_new)
	451
[101]	452	endif
[105]	453
[112]	454	tsurf_temp = tsurf_new
[98]	455	cal = 0.
	456	beta = 1.
	457	dif_grnd = 0.
[258]	458	alb_neig(:) = 0.
	459	agesno(:) = 0.
[98]	460
[147]	461	call calcul_fluxs( klon, knon, nisurf, dtime, &
[105]	462	& tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
[441]	463	& precip_rain, precip_snow, snow, qsurf, &
[90]	464	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
	465	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
[98]	466	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
[177]	467
[458]	468	fder_prev = fder
	469	fder = fder_prev + dflux_s + dflux_l
	470
	471	iloc = maxloc(fder(1:klon))
	472	if (check.and.fder(iloc(1))> 0.) then
	473	WRITE(,)'** Debug fder**'
	474	WRITE(,)'max fder(',iloc(1),') = ',fder(iloc(1))
	475	WRITE(,)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &
	476	& dflux_s(iloc(1)), dflux_l(iloc(1))
	477	endif
	478	!!$
	479	!!$ where(fder.gt.0.)
	480	!!$ fder = 0.
	481	!!$ endwhere
	482
[90]	483	!
[177]	484	! 2eme appel a interfoce pour le cumul des champs (en particulier
	485	! fluxsens et fluxlat calcules dans calcul_fluxs)
	486	!
	487	if (ocean == 'couple') then
	488
	489	cumul = .true.
	490
	491	call interfoce(itime, dtime, cumul, &
	492	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
	493	& ocean, npas, nexca, debut, lafin, &
	494	& swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
	495	& fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
[274]	496	& tsurf_new, alb_new, pctsrf_new)
[177]	497
	498	! else if (ocean == 'slab ') then
	499	! call interfoce(nisurf)
	500
	501	endif
	502
	503	!
[98]	504	! calcul albedo
	505	!
	506
[112]	507	if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then
	508	CALL alboc(FLOAT(jour),rlat,alb_eau)
	509	else ! cycle diurne
	510	CALL alboc_cd(rmu0,alb_eau)
	511	endif
	512	DO ii =1, knon
	513	alb_new(ii) = alb_eau(knindex(ii))
	514	enddo
[157]	515
[274]	516	z0_new = sqrt(rugos2 + rugoro2)
[281]	517	alblw(1:knon) = alb_new(1:knon)
	518
[98]	519	!
[90]	520	else if (nisurf == is_sic) then
	521
	522	if (check) write(,)'sea ice, nisurf = ',nisurf
	523
	524	!
	525	! Surface "glace de mer" appel a l'interface avec l'ocean
	526	!
	527	!
[105]	528	if (ocean == 'couple') then
[98]	529
[177]	530	cumul =.false.
	531
[458]	532	iloc = maxloc(fder(1:klon))
	533	if (check.and.fder(iloc(1))> 0.) then
	534	WRITE(,)'** Debug fder **'
	535	WRITE(,)'max fder(',iloc(1),') = ',fder(iloc(1))
	536	endif
	537	!!$
	538	!!$ where(fder.gt.0.)
	539	!!$ fder = 0.
	540	!!$ endwhere
	541
[177]	542	call interfoce(itime, dtime, cumul, &
[105]	543	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
[112]	544	& ocean, npas, nexca, debut, lafin, &
[177]	545	& swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
	546	& fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
[274]	547	& tsurf_new, alb_new, pctsrf_new)
[98]	548
[105]	549	tsurf_temp = tsurf_new
[177]	550	cal = 0.
[105]	551	dif_grnd = 0.
[140]	552	beta = 1.0
[105]	553
	554	! else if (ocean == 'slab ') then
	555	! call interfoce(nisurf)
[274]	556	ELSE
[235]	557	! ! lecture conditions limites
[274]	558	CALL interfoce(itime, dtime, jour, &
[235]	559	& klon, nisurf, knon, knindex, &
	560	& debut, &
	561	& tsurf_new, pctsrf_new)
[105]	562
[457]	563	!IM cf LF
[461]	564	DO ii = 1, knon
	565	IF (pctsrf_new(ii,nisurf) < EPSFRA) then
	566	snow(ii) = 0.0
	567	!IM cf LF/JLD tsurf(ii) = RTT - 1.8
	568	tsurf_new(ii) = RTT - 1.8
	569	IF (soil_model) tsoil(ii,:) = RTT -1.8
[457]	570	endif
	571	enddo
	572
[274]	573	CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
[235]	574
[274]	575	IF (soil_model) THEN
[457]	576	!IM cf LF/JLD CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)
	577	CALL soil(dtime, nisurf, knon,snow, tsurf_new, tsoil,soilcap, soilflux)
[274]	578	cal(1:knon) = RCPD / soilcap(1:knon)
	579	radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
	580	dif_grnd = 0.
	581	ELSE
	582	dif_grnd = 1.0 / tau_gl
	583	cal = RCPD * calice
	584	WHERE (snow > 0.0) cal = RCPD * calsno
[235]	585	ENDIF
[274]	586	tsurf_temp = tsurf
	587	beta = 1.0
	588	ENDIF
[105]	589
[274]	590	CALL calcul_fluxs( klon, knon, nisurf, dtime, &
[235]	591	& tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
[441]	592	& precip_rain, precip_snow, snow, qsurf, &
[235]	593	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
	594	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
	595	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
[90]	596
[457]	597	!IM cf JP 12.02.03
	598	! DO i = 1, knon
	599	! IF (pctsrf_new(i,nisurf) < EPSFRA) then
	600	! snow(i) = 0.0
	601	! tsurf_new(i) = RTT - 1.8
	602	! IF (soil_model) tsoil(i,:) = RTT -1.8
	603	! endif
	604	! enddo
	605
[274]	606	IF (ocean /= 'couple') THEN
	607	CALL fonte_neige( klon, knon, nisurf, dtime, &
[235]	608	& tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
	609	& precip_rain, precip_snow, snow, qsol, &
	610	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
	611	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
[457]	612	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
	613	!IM cf JLD
[522]	614	& fqcalving,ffonte, run_off_lic_0)
[274]	615
	616	! calcul albedo
	617
	618	CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
	619	WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
[295]	620	zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
[281]	621	alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
	622	& 0.6 * (1.0-zfra(1:knon))
[274]	623	!! alb_new(1 : knon) = 0.6
	624	ENDIF
	625
[458]	626	fder_prev = fder
	627	fder = fder_prev + dflux_s + dflux_l
[274]	628
[458]	629	iloc = maxloc(fder(1:klon))
	630	if (check.and.fder(iloc(1))> 0.) then
	631	WRITE(,)'** Debug fder **'
	632	WRITE(,)'max fder(',iloc(1),') = ',fder(iloc(1))
	633	WRITE(,)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &
	634	& dflux_s(iloc(1)), dflux_l(iloc(1))
	635	endif
	636	!!$ where(fder.gt.0.)
	637	!!$ fder = 0.
	638	!!$ endwhere
	639
[109]	640	!
[177]	641	! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
	642	!
	643	if (ocean == 'couple') then
	644
	645	cumul =.true.
	646
	647	call interfoce(itime, dtime, cumul, &
	648	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
	649	& ocean, npas, nexca, debut, lafin, &
	650	& swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
	651	& fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
[274]	652	& tsurf_new, alb_new, pctsrf_new)
[177]	653
	654	! else if (ocean == 'slab ') then
	655	! call interfoce(nisurf)
	656
	657	endif
	658
[258]	659
[489]	660	z0_new = 0.002
[281]	661	z0_new = SQRT(z0_new2+rugoro2)
	662	alblw(1:knon) = alb_new(1:knon)
[98]	663
[90]	664	else if (nisurf == is_lic) then
	665
	666	if (check) write(,)'glacier, nisurf = ',nisurf
	667
[522]	668	if (.not. allocated(run_off_lic)) then
	669	allocate(run_off_lic(knon), stat = error)
	670	if (error /= 0) then
	671	abort_message='Pb allocation run_off_lic'
	672	call abort_gcm(modname,abort_message,1)
	673	endif
	674	run_off_lic = 0.
	675	endif
[90]	676	!
	677	! Surface "glacier continentaux" appel a l'interface avec le sol
	678	!
	679	! call interfsol(nisurf)
[177]	680	IF (soil_model) THEN
[235]	681	CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil,soilcap, soilflux)
	682	cal(1:knon) = RCPD / soilcap(1:knon)
	683	radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
[177]	684	ELSE
	685	cal = RCPD * calice
	686	WHERE (snow > 0.0) cal = RCPD * calsno
	687	ENDIF
[98]	688	beta = 1.0
	689	dif_grnd = 0.0
	690
[147]	691	call calcul_fluxs( klon, knon, nisurf, dtime, &
[105]	692	& tsurf, p1lay, cal, beta, tq_cdrag, ps, &
[441]	693	& precip_rain, precip_snow, snow, qsurf, &
[90]	694	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
	695	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
[98]	696	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
[90]	697
[181]	698	call fonte_neige( klon, knon, nisurf, dtime, &
	699	& tsurf, p1lay, cal, beta, tq_cdrag, ps, &
	700	& precip_rain, precip_snow, snow, qsol, &
	701	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
	702	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
[457]	703	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
	704	!IM cf JLD
[522]	705	& fqcalving,ffonte, run_off_lic_0)
[181]	706
[522]	707	! passage du run-off des glaciers calcule dans fonte_neige au coupleur
	708	bidule=0.
	709	bidule(1:knon)= run_off_lic(1:knon)
	710	call gath2cpl(bidule, tmp_rlic, klon, knon,iim,jjm,knindex)
[109]	711	!
	712	! calcul albedo
	713	!
[258]	714	CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))
	715	WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
[295]	716	zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
[281]	717	alb_new(1 : knon) = alb_neig(1 : knon)*zfra(1:knon) + &
	718	& 0.6 * (1.0-zfra(1:knon))
[486]	719	!IM cf FH/GK alb_new(1 : knon) = 0.6
	720	! alb_new(1 : knon) = 0.82
	721	!IM cf JLD/ GK
	722	!IM: 211003 Ksta0.77 alb_new(1 : knon) = 0.77
	723	!IM: KstaTER0.8 & LMD_ARMIP5 alb_new(1 : knon) = 0.8
	724	!IM: KstaTER0.77 & LMD_ARMIP6
[513]	725	alb_new(1 : knon) = 0.77
[486]	726
[139]	727	!
[157]	728	! Rugosite
	729	!
[171]	730	z0_new = rugoro
[157]	731	!
[139]	732	! Remplissage des pourcentages de surface
	733	!
	734	pctsrf_new(:,nisurf) = pctsrf(:,nisurf)
[109]	735
[281]	736	alblw(1:knon) = alb_new(1:knon)
[90]	737	else
	738	write(,)'Index surface = ',nisurf
	739	abort_message = 'Index surface non valable'
	740	call abort_gcm(modname,abort_message,1)
	741	endif
	742
	743	END SUBROUTINE interfsurf_hq
	744
	745	!
	746	!#########################################################################
	747	!
	748	SUBROUTINE interfsurf_vent(nisurf, knon &
	749	& )
	750	!
	751	! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general
	752	! (sols continentaux, oceans, glaces) pour les tensions de vents.
	753	! En pratique l'interface se fait entre la couche limite du modele
	754	! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
	755	!
	756	!
	757	! L.Fairhead 02/2000
	758	!
	759	! input:
	760	! nisurf index de la surface a traiter (1 = sol continental)
	761	! knon nombre de points de la surface a traiter
	762
	763	! Parametres d'entree
	764	integer, intent(IN) :: nisurf
	765	integer, intent(IN) :: knon
	766
	767
	768	return
	769	END SUBROUTINE interfsurf_vent
	770	!
	771	!#########################################################################
	772	!
[205]	773	SUBROUTINE interfsol(itime, klon, dtime, date0, nisurf, knon, &
[177]	774	& knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
[90]	775	& debut, lafin, ok_veget, &
[441]	776	& plev, u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
[90]	777	& tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
	778	& precip_rain, precip_snow, lwdown, swnet, swdown, &
[105]	779	& tsurf, p1lay, ps, radsol, &
[90]	780	& evap, fluxsens, fluxlat, &
[281]	781	& tsol_rad, tsurf_new, alb_new, alblw, &
[365]	782	& emis_new, z0_new, dflux_l, dflux_s, qsurf)
[90]	783
[177]	784	USE intersurf
	785
[90]	786	! Cette routine sert d'interface entre le modele atmospherique et le
	787	! modele de sol continental. Appel a sechiba
	788	!
	789	! L. Fairhead 02/2000
	790	!
	791	! input:
	792	! itime numero du pas de temps
	793	! klon nombre total de points de grille
	794	! dtime pas de temps de la physique (en s)
	795	! nisurf index de la surface a traiter (1 = sol continental)
	796	! knon nombre de points de la surface a traiter
	797	! knindex index des points de la surface a traiter
	798	! rlon longitudes de la grille entiere
	799	! rlat latitudes de la grille entiere
[177]	800	! pctsrf tableau des fractions de surface de chaque maille
[90]	801	! debut logical: 1er appel a la physique (lire les restart)
	802	! lafin logical: dernier appel a la physique (ecrire les restart)
	803	! ok_veget logical: appel ou non au schema de surface continental
	804	! (si false calcul simplifie des fluxs sur les continents)
[441]	805	! plev hauteur de la premiere couche (Pa)
[90]	806	! u1_lay vitesse u 1ere couche
	807	! v1_lay vitesse v 1ere couche
	808	! temp_air temperature de l'air 1ere couche
	809	! spechum humidite specifique 1ere couche
[177]	810	! epot_air temp pot de l'air
[90]	811	! ccanopy concentration CO2 canopee
	812	! tq_cdrag cdrag
	813	! petAcoef coeff. A de la resolution de la CL pour t
	814	! peqAcoef coeff. A de la resolution de la CL pour q
	815	! petBcoef coeff. B de la resolution de la CL pour t
	816	! peqBcoef coeff. B de la resolution de la CL pour q
	817	! precip_rain precipitation liquide
	818	! precip_snow precipitation solide
[177]	819	! lwdown flux IR descendant a la surface
[90]	820	! swnet flux solaire net
	821	! swdown flux solaire entrant a la surface
	822	! tsurf temperature de surface
	823	! p1lay pression 1er niveau (milieu de couche)
	824	! ps pression au sol
	825	! radsol rayonnement net aus sol (LW + SW)
	826	!
	827	!
	828	! input/output
	829	! run_off ruissellement total
	830	!
	831	! output:
	832	! evap evaporation totale
	833	! fluxsens flux de chaleur sensible
	834	! fluxlat flux de chaleur latente
	835	! tsol_rad
	836	! tsurf_new temperature au sol
	837	! alb_new albedo
	838	! emis_new emissivite
	839	! z0_new surface roughness
[441]	840	! qsurf air moisture at surface
[90]	841
	842	! Parametres d'entree
	843	integer, intent(IN) :: itime
	844	integer, intent(IN) :: klon
	845	real, intent(IN) :: dtime
[205]	846	real, intent(IN) :: date0
[90]	847	integer, intent(IN) :: nisurf
	848	integer, intent(IN) :: knon
[177]	849	integer, intent(IN) :: iim, jjm
[147]	850	integer, dimension(klon), intent(IN) :: knindex
[90]	851	logical, intent(IN) :: debut, lafin, ok_veget
[177]	852	real, dimension(klon,nbsrf), intent(IN) :: pctsrf
[90]	853	real, dimension(klon), intent(IN) :: rlon, rlat
[177]	854	real, dimension(klon), intent(IN) :: cufi, cvfi
[441]	855	real, dimension(klon), intent(IN) :: plev
[147]	856	real, dimension(klon), intent(IN) :: u1_lay, v1_lay
	857	real, dimension(klon), intent(IN) :: temp_air, spechum
[177]	858	real, dimension(klon), intent(IN) :: epot_air, ccanopy
	859	real, dimension(klon), intent(INOUT) :: tq_cdrag
	860	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
[147]	861	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
	862	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
	863	real, dimension(klon), intent(IN) :: lwdown, swnet, swdown, ps
[441]	864	!IM cf. JP +++
	865	real, dimension(klon) :: swdown_vrai
	866	!IM cf. JP ---
[147]	867	real, dimension(klon), intent(IN) :: tsurf, p1lay
	868	real, dimension(klon), intent(IN) :: radsol
[90]	869	! Parametres de sortie
[365]	870	real, dimension(klon), intent(OUT):: evap, fluxsens, fluxlat, qsurf
[281]	871	real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new, alblw
[147]	872	real, dimension(klon), intent(OUT):: emis_new, z0_new
	873	real, dimension(klon), intent(OUT):: dflux_s, dflux_l
[90]	874
	875	! Local
	876	!
[281]	877	integer :: ii, ij, jj, igrid, ireal, i, index, iglob
[90]	878	integer :: error
	879	character (len = 20) :: modname = 'interfsol'
	880	character (len = 80) :: abort_message
[295]	881	logical,save :: check = .FALSE.
[147]	882	real, dimension(klon) :: cal, beta, dif_grnd, capsol
[90]	883	! type de couplage dans sechiba
	884	! character (len=10) :: coupling = 'implicit'
	885	! drapeaux controlant les appels dans SECHIBA
	886	! type(control_type), save :: control_in
[441]	887	! Preserved albedo
	888	!IM cf. JP +++
	889	real, allocatable, dimension(:), save :: albedo_keep, zlev
	890	!IM cf. JP ---
[90]	891	! coordonnees geographiques
	892	real, allocatable, dimension(:,:), save :: lalo
	893	! pts voisins
	894	integer,allocatable, dimension(:,:), save :: neighbours
[177]	895	! fractions continents
	896	real,allocatable, dimension(:), save :: contfrac
[90]	897	! resolution de la grille
	898	real, allocatable, dimension (:,:), save :: resolution
[177]	899	! correspondance point n -> indices (i,j)
	900	integer, allocatable, dimension(:,:), save :: correspond
	901	! offset pour calculer les point voisins
	902	integer, dimension(8,3), save :: off_ini
	903	integer, dimension(8), save :: offset
[90]	904	! Identifieurs des fichiers restart et histoire
	905	integer, save :: rest_id, hist_id
	906	integer, save :: rest_id_stom, hist_id_stom
[177]	907	!
	908	real, allocatable, dimension (:,:), save :: lon_scat, lat_scat
[90]	909
[223]	910	logical, save :: lrestart_read = .true. , lrestart_write = .false.
[177]	911
[365]	912	real, dimension(klon):: snow
[177]	913	real, dimension(knon,2) :: albedo_out
	914	! Pb de nomenclature
	915	real, dimension(klon) :: petA_orc, peqA_orc
	916	real, dimension(klon) :: petB_orc, peqB_orc
[281]	917	! Pb de correspondances de grilles
	918	integer, dimension(:), save, allocatable :: ig, jg
	919	integer :: indi, indj
	920	integer, dimension(klon) :: ktindex
[290]	921	REAL, dimension(klon) :: bidule
[281]	922	! Essai cdrag
	923	real, dimension(klon) :: cdrag
[98]	924
[397]	925	#include "temps.inc"
[441]	926	#include "YOMCST.inc"
[364]	927
[90]	928	if (check) write(,)'Entree ', modname
	929	if (check) write(,)'ok_veget = ',ok_veget
	930
[281]	931	ktindex(:) = knindex(:) + iim - 1
	932
[90]	933	! initialisation
[177]	934	if (debut) then
[90]	935
[441]	936	IF ( .NOT. allocated(albedo_keep)) THEN
	937	ALLOCATE(albedo_keep(klon))
	938	ALLOCATE(zlev(klon))
	939	ENDIF
[281]	940	! Pb de correspondances de grilles
	941	allocate(ig(klon))
	942	allocate(jg(klon))
	943	ig(1) = 1
	944	jg(1) = 1
	945	indi = 0
	946	indj = 2
	947	do igrid = 2, klon - 1
	948	indi = indi + 1
	949	if ( indi > iim) then
	950	indi = 1
	951	indj = indj + 1
	952	endif
	953	ig(igrid) = indi
	954	jg(igrid) = indj
	955	enddo
	956	ig(klon) = 1
	957	jg(klon) = jjm + 1
[177]	958	!
	959	! Initialisation des offset
	960	!
	961	! offset bord ouest
	962	off_ini(1,1) = - iim ; off_ini(2,1) = - iim + 1; off_ini(3,1) = 1
	963	off_ini(4,1) = iim + 1; off_ini(5,1) = iim ; off_ini(6,1) = 2 * iim - 1
	964	off_ini(7,1) = iim -1 ; off_ini(8,1) = - 1
	965	! offset point normal
	966	off_ini(1,2) = - iim ; off_ini(2,2) = - iim + 1; off_ini(3,2) = 1
	967	off_ini(4,2) = iim + 1; off_ini(5,2) = iim ; off_ini(6,2) = iim - 1
	968	off_ini(7,2) = -1 ; off_ini(8,2) = - iim - 1
	969	! offset bord est
	970	off_ini(1,3) = - iim; off_ini(2,3) = - 2 * iim + 1; off_ini(3,3) = - iim + 1
	971	off_ini(4,3) = 1 ; off_ini(5,3) = iim ; off_ini(6,3) = iim - 1
	972	off_ini(7,3) = -1 ; off_ini(8,3) = - iim - 1
	973	!
	974	! Initialisation des correspondances point -> indices i,j
	975	!
	976	if (( .not. allocated(correspond))) then
	977	allocate(correspond(iim,jjm+1), stat = error)
	978	if (error /= 0) then
	979	abort_message='Pb allocation correspond'
	980	call abort_gcm(modname,abort_message,1)
	981	endif
	982	endif
	983	!
	984	! Attention aux poles
	985	!
	986	do igrid = 1, knon
[281]	987	index = ktindex(igrid)
	988	jj = int((index - 1)/iim) + 1
	989	ij = index - (jj - 1) * iim
[177]	990	correspond(ij,jj) = igrid
	991	enddo
[281]	992
[177]	993	! Allouer et initialiser le tableau de coordonnees du sol
	994	!
	995	if ((.not. allocated(lalo))) then
	996	allocate(lalo(knon,2), stat = error)
	997	if (error /= 0) then
	998	abort_message='Pb allocation lalo'
	999	call abort_gcm(modname,abort_message,1)
	1000	endif
	1001	endif
	1002	if ((.not. allocated(lon_scat))) then
[201]	1003	allocate(lon_scat(iim,jjm+1), stat = error)
[177]	1004	if (error /= 0) then
	1005	abort_message='Pb allocation lon_scat'
	1006	call abort_gcm(modname,abort_message,1)
	1007	endif
	1008	endif
	1009	if ((.not. allocated(lat_scat))) then
[201]	1010	allocate(lat_scat(iim,jjm+1), stat = error)
[177]	1011	if (error /= 0) then
	1012	abort_message='Pb allocation lat_scat'
	1013	call abort_gcm(modname,abort_message,1)
	1014	endif
	1015	endif
	1016	lon_scat = 0.
	1017	lat_scat = 0.
	1018	do igrid = 1, knon
	1019	index = knindex(igrid)
	1020	lalo(igrid,2) = rlon(index)
	1021	lalo(igrid,1) = rlat(index)
	1022	ij = index - int((index-1)/iim)*iim - 1
	1023	jj = 2 + int((index-1)/iim)
	1024	if (mod(index,iim) == 1 ) then
	1025	jj = 1 + int((index-1)/iim)
	1026	ij = iim
	1027	endif
[281]	1028	! lon_scat(ij,jj) = rlon(index)
	1029	! lat_scat(ij,jj) = rlat(index)
[177]	1030	enddo
	1031	index = 1
	1032	do jj = 2, jjm
	1033	do ij = 1, iim
	1034	index = index + 1
	1035	lon_scat(ij,jj) = rlon(index)
	1036	lat_scat(ij,jj) = rlat(index)
	1037	enddo
	1038	enddo
	1039	lon_scat(:,1) = lon_scat(:,2)
	1040	lat_scat(:,1) = rlat(1)
[201]	1041	lon_scat(:,jjm+1) = lon_scat(:,2)
	1042	lat_scat(:,jjm+1) = rlat(klon)
[281]	1043	! Pb de correspondances de grilles!
	1044	! do igrid = 1, knon
	1045	! index = ktindex(igrid)
	1046	! ij = ig(index)
	1047	! jj = jg(index)
	1048	! lon_scat(ij,jj) = rlon(index)
	1049	! lat_scat(ij,jj) = rlat(index)
	1050	! enddo
[90]	1051
[177]	1052	!
	1053	! Allouer et initialiser le tableau des voisins et des fraction de continents
	1054	!
	1055	if ( (.not.allocated(neighbours))) THEN
	1056	allocate(neighbours(knon,8), stat = error)
	1057	if (error /= 0) then
	1058	abort_message='Pb allocation neighbours'
	1059	call abort_gcm(modname,abort_message,1)
	1060	endif
	1061	endif
[274]	1062	neighbours = -1.
[177]	1063	if (( .not. allocated(contfrac))) then
	1064	allocate(contfrac(knon), stat = error)
	1065	if (error /= 0) then
	1066	abort_message='Pb allocation contfrac'
	1067	call abort_gcm(modname,abort_message,1)
	1068	endif
	1069	endif
[90]	1070
[177]	1071	do igrid = 1, knon
	1072	ireal = knindex(igrid)
	1073	contfrac(igrid) = pctsrf(ireal,is_ter)
[281]	1074	enddo
	1075
	1076	do igrid = 1, knon
	1077	iglob = ktindex(igrid)
	1078	if (mod(iglob, iim) == 1) then
[177]	1079	offset = off_ini(:,1)
[281]	1080	else if(mod(iglob, iim) == 0) then
[177]	1081	offset = off_ini(:,3)
	1082	else
	1083	offset = off_ini(:,2)
	1084	endif
	1085	do i = 1, 8
[281]	1086	index = iglob + offset(i)
	1087	ireal = (min(max(1, index - iim + 1), klon))
	1088	if (pctsrf(ireal, is_ter) > EPSFRA) then
	1089	jj = int((index - 1)/iim) + 1
	1090	ij = index - (jj - 1) * iim
[177]	1091	neighbours(igrid, i) = correspond(ij, jj)
	1092	endif
	1093	enddo
	1094	enddo
	1095
	1096	!
	1097	! Allocation et calcul resolutions
	1098	IF ( (.NOT.ALLOCATED(resolution))) THEN
	1099	ALLOCATE(resolution(knon,2), stat = error)
	1100	if (error /= 0) then
	1101	abort_message='Pb allocation resolution'
	1102	call abort_gcm(modname,abort_message,1)
	1103	endif
	1104	ENDIF
	1105	do igrid = 1, knon
	1106	ij = knindex(igrid)
	1107	resolution(igrid,1) = cufi(ij)
	1108	resolution(igrid,2) = cvfi(ij)
	1109	enddo
	1110
	1111	endif ! (fin debut)
	1112
[90]	1113	!
	1114	! Appel a la routine sols continentaux
	1115	!
[223]	1116	if (lafin) lrestart_write = .true.
	1117	if (check) write(,)'lafin ',lafin,lrestart_write
[90]	1118
[177]	1119	petA_orc = petBcoef * dtime
	1120	petB_orc = petAcoef
	1121	peqA_orc = peqBcoef * dtime
	1122	peqB_orc = peqAcoef
[90]	1123
[281]	1124	cdrag = 0.
	1125	cdrag(1:knon) = tq_cdrag(1:knon)
	1126
[441]	1127	!IM cf. JP +++
	1128	! zlev(1:knon) = (100.plev(1:knon))/((ps(1:knon)/287.05temp_air(1:knon))*9.80665)
	1129	zlev(1:knon) = (100.plev(1:knon))/((ps(1:knon)/RDtemp_air(1:knon))*RG)
	1130	!IM cf. JP ---
	1131
[467]	1132
	1133	! PF et PASB
	1134	! where(cdrag > 0.01)
	1135	! cdrag = 0.01
	1136	! endwhere
[281]	1137	! write(,)'Cdrag = ',minval(cdrag),maxval(cdrag)
	1138
[201]	1139	!
	1140	! Init Orchidee
	1141	!
	1142	if (debut) then
[364]	1143	call intersurf_main (itime+itau_phy-1, iim, jjm+1, knon, ktindex, dtime, &
[177]	1144	& lrestart_read, lrestart_write, lalo, &
	1145	& contfrac, neighbours, resolution, date0, &
	1146	& zlev, u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
[281]	1147	& cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
[360]	1148	& precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
[177]	1149	& evap, fluxsens, fluxlat, coastalflow, riverflow, &
	1150	& tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
	1151	& lon_scat, lat_scat)
[441]	1152
	1153	!IM cf. JP +++
[467]	1154	albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
[441]	1155	!IM cf. JP ---
	1156
[201]	1157	endif
[90]	1158
[441]	1159	!IM cf. JP +++
[467]	1160	swdown_vrai(1:knon) = swnet(1:knon)/(1. - albedo_keep(1:knon))
[441]	1161	!IM cf. JP ---
	1162
[364]	1163	call intersurf_main (itime+itau_phy, iim, jjm+1, knon, ktindex, dtime, &
[201]	1164	& lrestart_read, lrestart_write, lalo, &
	1165	& contfrac, neighbours, resolution, date0, &
	1166	& zlev, u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
[281]	1167	& cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
[441]	1168	!IM cf. JP +++
	1169	& precip_rain, precip_snow, lwdown, swnet, swdown_vrai, ps, &
	1170	!IM cf. JP ---
[201]	1171	& evap, fluxsens, fluxlat, coastalflow, riverflow, &
	1172	& tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
	1173	& lon_scat, lat_scat)
	1174
[441]	1175	!IM cf. JP +++
[486]	1176	!IM BUG BUG BUG albedo_keep(:) = (albedo_out(:,1)+albedo_out(:,2))/2.
	1177	albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
[441]	1178	!IM cf. JP ---
	1179
[290]	1180	bidule=0.
	1181	bidule(1:knon)=riverflow(1:knon)
[398]	1182	call gath2cpl(bidule, tmp_rriv, klon, knon,iim,jjm,knindex)
[290]	1183	bidule=0.
	1184	bidule(1:knon)=coastalflow(1:knon)
[398]	1185	call gath2cpl(bidule, tmp_rcoa, klon, knon,iim,jjm,knindex)
[290]	1186	alb_new(1:knon) = albedo_out(1:knon,1)
	1187	alblw(1:knon) = albedo_out(1:knon,2)
[177]	1188
[290]	1189
[275]	1190	! Convention orchidee: positif vers le haut
[290]	1191	fluxsens(1:knon) = -1. * fluxsens(1:knon)
	1192	fluxlat(1:knon) = -1. * fluxlat(1:knon)
	1193
[281]	1194	! evap = -1. * evap
[223]	1195
	1196	if (debut) lrestart_read = .false.
	1197
[90]	1198	END SUBROUTINE interfsol
	1199	!
	1200	!#########################################################################
	1201	!
[177]	1202	SUBROUTINE interfoce_cpl(itime, dtime, cumul, &
[98]	1203	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
[112]	1204	& ocean, npas, nexca, debut, lafin, &
[98]	1205	& swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
[177]	1206	& fluxlat, fluxsens, fder, albsol, taux, tauy, zmasq, &
[274]	1207	& tsurf_new, alb_new, pctsrf_new)
[90]	1208
	1209	! Cette routine sert d'interface entre le modele atmospherique et un
[101]	1210	! coupleur avec un modele d'ocean 'complet' derriere
[90]	1211	!
[105]	1212	! Le modele de glace qu'il est prevu d'utiliser etant couple directement a
	1213	! l'ocean presentement, on va passer deux fois dans cette routine par pas de
	1214	! temps physique, une fois avec les points oceans et l'autre avec les points
	1215	! glace. A chaque pas de temps de couplage, la lecture des champs provenant
	1216	! du coupleur se fera "dans" l'ocean et l'ecriture des champs a envoyer
	1217	! au coupleur "dans" la glace. Il faut donc des tableaux de travail "tampons"
	1218	! dimensionnes sur toute la grille qui remplissent les champs sur les
	1219	! domaines ocean/glace quand il le faut. Il est aussi necessaire que l'index
	1220	! ocean soit traiter avant l'index glace (sinon tout intervertir)
	1221	!
	1222	!
[90]	1223	! L. Fairhead 02/2000
	1224	!
	1225	! input:
[98]	1226	! itime numero du pas de temps
	1227	! iim, jjm nbres de pts de grille
[313]	1228	! dtime pas de temps de la physique
[98]	1229	! klon nombre total de points de grille
[90]	1230	! nisurf index de la surface a traiter (1 = sol continental)
[98]	1231	! pctsrf tableau des fractions de surface de chaque maille
	1232	! knon nombre de points de la surface a traiter
	1233	! knindex index des points de la surface a traiter
	1234	! rlon longitudes
	1235	! rlat latitudes
	1236	! debut logical: 1er appel a la physique
	1237	! lafin logical: dernier appel a la physique
[90]	1238	! ocean type d'ocean
[98]	1239	! nexca frequence de couplage
	1240	! swdown flux solaire entrant a la surface
[177]	1241	! lwdown flux IR net a la surface
[98]	1242	! precip_rain precipitation liquide
	1243	! precip_snow precipitation solide
	1244	! evap evaporation
	1245	! tsurf temperature de surface
	1246	! fder derivee dF/dT
	1247	! albsol albedo du sol (coherent avec swdown)
	1248	! taux tension de vent en x
	1249	! tauy tension de vent en y
	1250	! nexca frequence de couplage
[105]	1251	! zmasq masque terre/ocean
[90]	1252	!
[98]	1253	!
[90]	1254	! output:
[98]	1255	! tsurf_new temperature au sol
	1256	! alb_new albedo
	1257	! pctsrf_new nouvelle repartition des surfaces
	1258	! alb_ice albedo de la glace
[90]	1259	!
	1260
[98]	1261
[90]	1262	! Parametres d'entree
[98]	1263	integer, intent(IN) :: itime
	1264	integer, intent(IN) :: iim, jjm
	1265	real, intent(IN) :: dtime
	1266	integer, intent(IN) :: klon
[90]	1267	integer, intent(IN) :: nisurf
[98]	1268	integer, intent(IN) :: knon
	1269	real, dimension(klon,nbsrf), intent(IN) :: pctsrf
[147]	1270	integer, dimension(klon), intent(in) :: knindex
[98]	1271	logical, intent(IN) :: debut, lafin
	1272	real, dimension(klon), intent(IN) :: rlon, rlat
[90]	1273	character (len = 6) :: ocean
[147]	1274	real, dimension(klon), intent(IN) :: lwdown, swdown
	1275	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
	1276	real, dimension(klon), intent(IN) :: tsurf, fder, albsol, taux, tauy
[140]	1277	INTEGER :: nexca, npas, kstep
[105]	1278	real, dimension(klon), intent(IN) :: zmasq
[177]	1279	real, dimension(klon), intent(IN) :: fluxlat, fluxsens
	1280	logical, intent(IN) :: cumul
[147]	1281	real, dimension(klon), intent(INOUT) :: evap
[98]	1282
[90]	1283	! Parametres de sortie
[274]	1284	real, dimension(klon), intent(OUT):: tsurf_new, alb_new
[98]	1285	real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
[90]	1286
	1287	! Variables locales
[140]	1288	integer :: j, error, sum_error, ig, cpl_index,i
[98]	1289	character (len = 20) :: modname = 'interfoce_cpl'
	1290	character (len = 80) :: abort_message
[295]	1291	logical,save :: check = .FALSE.
[98]	1292	! variables pour moyenner les variables de couplage
[105]	1293	real, allocatable, dimension(:,:),save :: cpl_sols, cpl_nsol, cpl_rain
	1294	real, allocatable, dimension(:,:),save :: cpl_snow, cpl_evap, cpl_tsol
	1295	real, allocatable, dimension(:,:),save :: cpl_fder, cpl_albe, cpl_taux
[290]	1296	real, allocatable, dimension(:,:),save :: cpl_tauy
[522]	1297	REAL, ALLOCATABLE, DIMENSION(:,:),SAVE :: cpl_rriv, cpl_rcoa, cpl_rlic
[290]	1298	!!$
[105]	1299	! variables tampons avant le passage au coupleur
	1300	real, allocatable, dimension(:,:,:),save :: tmp_sols, tmp_nsol, tmp_rain
	1301	real, allocatable, dimension(:,:,:),save :: tmp_snow, tmp_evap, tmp_tsol
	1302	real, allocatable, dimension(:,:,:),save :: tmp_fder, tmp_albe, tmp_taux
[398]	1303	!!$ real, allocatable, dimension(:,:,:),save :: tmp_tauy, tmp_rriv, tmp_rcoa
	1304	REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: tmp_tauy
[98]	1305	! variables a passer au coupleur
[105]	1306	real, dimension(iim, jjm+1) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice
	1307	real, dimension(iim, jjm+1) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
[177]	1308	REAL, DIMENSION(iim, jjm+1) :: wri_evap_sea, wri_rcoa, wri_rriv
	1309	REAL, DIMENSION(iim, jjm+1) :: wri_rain, wri_snow, wri_taux, wri_tauy
[394]	1310	REAL, DIMENSION(iim, jjm+1) :: wri_calv
[177]	1311	REAL, DIMENSION(iim, jjm+1) :: wri_tauxx, wri_tauyy, wri_tauzz
	1312	REAL, DIMENSION(iim, jjm+1) :: tmp_lon, tmp_lat
[98]	1313	! variables relues par le coupleur
[105]	1314	! read_sic = fraction de glace
	1315	! read_sit = temperature de glace
	1316	real, allocatable, dimension(:,:),save :: read_sst, read_sic, read_sit
	1317	real, allocatable, dimension(:,:),save :: read_alb_sic
[98]	1318	! variable tampon
[140]	1319	real, dimension(klon) :: tamp_sic
	1320	! sauvegarde des fractions de surface d'un pas de temps a l'autre apres
	1321	! l'avoir lu
	1322	real, allocatable,dimension(:,:),save :: pctsrf_sav
[105]	1323	real, dimension(iim, jjm+1, 2) :: tamp_srf
	1324	integer, allocatable, dimension(:), save :: tamp_ind
	1325	real, allocatable, dimension(:,:),save :: tamp_zmasq
	1326	real, dimension(iim, jjm+1) :: deno
[137]	1327	integer :: idtime
[179]	1328	integer, allocatable,dimension(:),save :: unity
[131]	1329	!
	1330	logical, save :: first_appel = .true.
[179]	1331	logical,save :: print
	1332	!maf
	1333	! variables pour avoir une sortie IOIPSL des champs echanges
	1334	CHARACTER*80,SAVE :: clintocplnam, clfromcplnam
	1335	INTEGER, SAVE :: jf,nhoridct,nidct
	1336	INTEGER, SAVE :: nhoridcs,nidcs
	1337	INTEGER :: ndexct(iim(jjm+1)),ndexcs(iim(jjm+1))
	1338	REAL :: zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
[353]	1339	integer :: idayref, itau_w
[397]	1340	#include "param_cou.h"
	1341	#include "inc_cpl.h"
	1342	#include "temps.inc"
[98]	1343	!
[90]	1344	! Initialisation
[98]	1345	!
[131]	1346	if (check) write(,)'Entree ',modname,'nisurf = ',nisurf
	1347
	1348	if (first_appel) then
[140]	1349	error = 0
	1350	allocate(unity(klon), stat = error)
	1351	if ( error /=0) then
	1352	abort_message='Pb allocation variable unity'
	1353	call abort_gcm(modname,abort_message,1)
	1354	endif
[143]	1355	allocate(pctsrf_sav(klon,nbsrf), stat = error)
[140]	1356	if ( error /=0) then
	1357	abort_message='Pb allocation variable pctsrf_sav'
	1358	call abort_gcm(modname,abort_message,1)
	1359	endif
[143]	1360	pctsrf_sav = 0.
[140]	1361
	1362	do ig = 1, klon
	1363	unity(ig) = ig
	1364	enddo
[98]	1365	sum_error = 0
[147]	1366	allocate(cpl_sols(klon,2), stat = error); sum_error = sum_error + error
	1367	allocate(cpl_nsol(klon,2), stat = error); sum_error = sum_error + error
	1368	allocate(cpl_rain(klon,2), stat = error); sum_error = sum_error + error
	1369	allocate(cpl_snow(klon,2), stat = error); sum_error = sum_error + error
	1370	allocate(cpl_evap(klon,2), stat = error); sum_error = sum_error + error
	1371	allocate(cpl_tsol(klon,2), stat = error); sum_error = sum_error + error
	1372	allocate(cpl_fder(klon,2), stat = error); sum_error = sum_error + error
	1373	allocate(cpl_albe(klon,2), stat = error); sum_error = sum_error + error
	1374	allocate(cpl_taux(klon,2), stat = error); sum_error = sum_error + error
	1375	allocate(cpl_tauy(klon,2), stat = error); sum_error = sum_error + error
[290]	1376	ALLOCATE(cpl_rriv(iim,jjm+1), stat=error); sum_error = sum_error + error
	1377	ALLOCATE(cpl_rcoa(iim,jjm+1), stat=error); sum_error = sum_error + error
[522]	1378	ALLOCATE(cpl_rlic(iim,jjm+1), stat=error); sum_error = sum_error + error
[290]	1379	!!
[105]	1380	allocate(read_sst(iim, jjm+1), stat = error); sum_error = sum_error + error
	1381	allocate(read_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
	1382	allocate(read_sit(iim, jjm+1), stat = error); sum_error = sum_error + error
	1383	allocate(read_alb_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
	1384
[98]	1385	if (sum_error /= 0) then
	1386	abort_message='Pb allocation variables couplees'
	1387	call abort_gcm(modname,abort_message,1)
	1388	endif
[105]	1389	cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
	1390	cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
[522]	1391	cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.; cpl_rlic = 0.
[105]	1392
	1393	sum_error = 0
	1394	allocate(tamp_ind(klon), stat = error); sum_error = sum_error + error
	1395	allocate(tamp_zmasq(iim, jjm+1), stat = error); sum_error = sum_error + error
	1396	do ig = 1, klon
	1397	tamp_ind(ig) = ig
	1398	enddo
	1399	call gath2cpl(zmasq, tamp_zmasq, klon, klon, iim, jjm, tamp_ind)
[98]	1400	!
	1401	! initialisation couplage
	1402	!
[137]	1403	idtime = int(dtime)
	1404	call inicma(npas , nexca, idtime,(jjm+1)*iim)
[98]	1405
[179]	1406	!
	1407	! initialisation sorties netcdf
	1408	!
[353]	1409	idayref = day_ini
	1410	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
[179]	1411	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
	1412	DO i = 1, iim
	1413	zx_lon(i,1) = rlon(i+1)
	1414	zx_lon(i,jjm+1) = rlon(i+1)
	1415	ENDDO
	1416	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
	1417	clintocplnam="cpl_atm_tauflx"
[353]	1418	CALL histbeg(clintocplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
	1419	& itau_phy,zjulian,dtime,nhoridct,nidct)
[179]	1420	! no vertical axis
	1421	CALL histdef(nidct, 'tauxe','tauxe', &
	1422	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
	1423	CALL histdef(nidct, 'tauyn','tauyn', &
	1424	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
	1425	CALL histdef(nidct, 'tmp_lon','tmp_lon', &
	1426	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
	1427	CALL histdef(nidct, 'tmp_lat','tmp_lat', &
	1428	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
	1429	DO jf=1,jpflda2o1 + jpflda2o2
	1430	CALL histdef(nidct, cl_writ(jf),cl_writ(jf), &
	1431	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
	1432	END DO
	1433	CALL histend(nidct)
	1434	CALL histsync(nidct)
	1435
	1436	clfromcplnam="cpl_atm_sst"
[353]	1437	CALL histbeg(clfromcplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
[179]	1438	& 0,zjulian,dtime,nhoridcs,nidcs)
	1439	! no vertical axis
	1440	DO jf=1,jpfldo2a
	1441	CALL histdef(nidcs, cl_read(jf),cl_read(jf), &
	1442	& "-",iim, jjm+1, nhoridcs, 1, 1, 1, -99, 32, "inst", dtime,dtime)
	1443	END DO
	1444	CALL histend(nidcs)
	1445	CALL histsync(nidcs)
	1446
[313]	1447	! pour simuler la fonte des glaciers antarctiques
	1448	!
	1449	surf_maille = (4. * rpi * ra*2) / (iim (jjm +1))
	1450	ALLOCATE(coeff_iceberg(iim,jjm+1), stat=error)
	1451	if (error /= 0) then
	1452	abort_message='Pb allocation variable coeff_iceberg'
	1453	call abort_gcm(modname,abort_message,1)
	1454	endif
	1455	open (12,file='flux_iceberg',form='formatted',status='old')
	1456	read (12,*) coeff_iceberg
	1457	close (12)
	1458	num_antarctic = max(1, count(coeff_iceberg > 0))
	1459
[131]	1460	first_appel = .false.
	1461	endif ! fin if (first_appel)
[98]	1462
[157]	1463	! Initialisations
[90]	1464
[112]	1465	! calcul des fluxs a passer
[90]	1466
[140]	1467	cpl_index = 1
	1468	if (nisurf == is_sic) cpl_index = 2
[177]	1469	if (cumul) then
[179]	1470	if (check) write(,) modname, 'cumul des champs'
[177]	1471	do ig = 1, knon
	1472	cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) &
	1473	& + swdown(ig) / FLOAT(nexca)
	1474	cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) &
	1475	& + (lwdown(ig) + fluxlat(ig) +fluxsens(ig))&
	1476	& / FLOAT(nexca)
	1477	cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) &
	1478	& + precip_rain(ig) / FLOAT(nexca)
	1479	cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) &
	1480	& + precip_snow(ig) / FLOAT(nexca)
	1481	cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) &
	1482	& + evap(ig) / FLOAT(nexca)
	1483	cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) &
	1484	& + tsurf(ig) / FLOAT(nexca)
	1485	cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) &
	1486	& + fder(ig) / FLOAT(nexca)
	1487	cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) &
	1488	& + albsol(ig) / FLOAT(nexca)
	1489	cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) &
	1490	& + taux(ig) / FLOAT(nexca)
	1491	cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) &
	1492	& + tauy(ig) / FLOAT(nexca)
	1493	enddo
[398]	1494	IF (cpl_index .EQ. 1) THEN
	1495	cpl_rriv(:,:) = cpl_rriv(:,:) + tmp_rriv(:,:) / FLOAT(nexca)
	1496	cpl_rcoa(:,:) = cpl_rcoa(:,:) + tmp_rcoa(:,:) / FLOAT(nexca)
[522]	1497	cpl_rlic(:,:) = cpl_rlic(:,:) + tmp_rlic(:,:) / FLOAT(nexca)
[398]	1498	ENDIF
[177]	1499	endif
[98]	1500
[138]	1501	if (mod(itime, nexca) == 1) then
[98]	1502	!
[179]	1503	! Demande des champs au coupleur
[98]	1504	!
[105]	1505	! Si le domaine considere est l'ocean, on lit les champs venant du coupleur
	1506	!
[177]	1507	if (nisurf == is_oce .and. .not. cumul) then
[147]	1508	if (check) write(,)'rentree fromcpl, itime-1 = ',itime-1
[138]	1509	call fromcpl(itime-1,(jjm+1)*iim, &
[105]	1510	& read_sst, read_sic, read_sit, read_alb_sic)
[179]	1511	!
	1512	! sorties NETCDF des champs recus
	1513	!
	1514	ndexcs(:)=0
[353]	1515	itau_w = itau_phy + itime
	1516	CALL histwrite(nidcs,cl_read(1),itau_w,read_sst,iim*(jjm+1),ndexcs)
	1517	CALL histwrite(nidcs,cl_read(2),itau_w,read_sic,iim*(jjm+1),ndexcs)
	1518	CALL histwrite(nidcs,cl_read(3),itau_w,read_alb_sic,iim*(jjm+1),ndexcs)
	1519	CALL histwrite(nidcs,cl_read(4),itau_w,read_sit,iim*(jjm+1),ndexcs)
[179]	1520	CALL histsync(nidcs)
	1521	! pas utile IF (npas-itime.LT.nexca )CALL histclo(nidcs)
	1522
[105]	1523	do j = 1, jjm + 1
	1524	do ig = 1, iim
	1525	if (abs(1. - read_sic(ig,j)) < 0.00001) then
	1526	read_sst(ig,j) = RTT - 1.8
	1527	read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
	1528	read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
	1529	else if (abs(read_sic(ig,j)) < 0.00001) then
	1530	read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
	1531	read_sit(ig,j) = read_sst(ig,j)
	1532	read_alb_sic(ig,j) = 0.6
	1533	else
	1534	read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
	1535	read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
	1536	read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
	1537	endif
	1538	enddo
	1539	enddo
[140]	1540	!
	1541	! transformer read_sic en pctsrf_sav
	1542	!
[157]	1543	call cpl2gath(read_sic, tamp_sic , klon, klon,iim,jjm, unity)
	1544	do ig = 1, klon
	1545	IF (pctsrf(ig,is_oce) > epsfra .OR. &
[140]	1546	& pctsrf(ig,is_sic) > epsfra) THEN
[143]	1547	pctsrf_sav(ig,is_sic) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
	1548	& * tamp_sic(ig)
	1549	pctsrf_sav(ig,is_oce) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
	1550	& - pctsrf_sav(ig,is_sic)
	1551	endif
	1552	enddo
[157]	1553	!
	1554	! Pour rattraper des erreurs d'arrondis
	1555	!
[179]	1556	where (abs(pctsrf_sav(:,is_sic)) .le. 2.*epsilon(pctsrf_sav(1,is_sic)))
[157]	1557	pctsrf_sav(:,is_sic) = 0.
	1558	pctsrf_sav(:,is_oce) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
	1559	endwhere
[179]	1560	where (abs(pctsrf_sav(:,is_oce)) .le. 2.*epsilon(pctsrf_sav(1,is_oce)))
[157]	1561	pctsrf_sav(:,is_sic) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
	1562	pctsrf_sav(:,is_oce) = 0.
	1563	endwhere
	1564	if (minval(pctsrf_sav(:,is_oce)) < 0.) then
[143]	1565	write(,)'Pb fraction ocean inferieure a 0'
[157]	1566	write(,)'au point ',minloc(pctsrf_sav(:,is_oce))
	1567	write(,)'valeur = ',minval(pctsrf_sav(:,is_oce))
[143]	1568	abort_message = 'voir ci-dessus'
	1569	call abort_gcm(modname,abort_message,1)
	1570	endif
[157]	1571	if (minval(pctsrf_sav(:,is_sic)) < 0.) then
[143]	1572	write(,)'Pb fraction glace inferieure a 0'
[157]	1573	write(,)'au point ',minloc(pctsrf_sav(:,is_sic))
	1574	write(,)'valeur = ',minval(pctsrf_sav(:,is_sic))
[143]	1575	abort_message = 'voir ci-dessus'
	1576	call abort_gcm(modname,abort_message,1)
	1577	endif
[105]	1578	endif
[138]	1579	endif ! fin mod(itime, nexca) == 1
	1580
	1581	if (mod(itime, nexca) == 0) then
[105]	1582	!
[138]	1583	! allocation memoire
[180]	1584	if (nisurf == is_oce .and. (.not. cumul) ) then
[147]	1585	sum_error = 0
	1586	allocate(tmp_sols(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1587	allocate(tmp_nsol(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1588	allocate(tmp_rain(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1589	allocate(tmp_snow(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1590	allocate(tmp_evap(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1591	allocate(tmp_tsol(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1592	allocate(tmp_fder(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1593	allocate(tmp_albe(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1594	allocate(tmp_taux(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1595	allocate(tmp_tauy(iim,jjm+1,2), stat=error); sum_error = sum_error + error
[398]	1596	!!$ allocate(tmp_rriv(iim,jjm+1,2), stat=error); sum_error = sum_error + error
	1597	!!$ allocate(tmp_rcoa(iim,jjm+1,2), stat=error); sum_error = sum_error + error
[147]	1598	if (sum_error /= 0) then
	1599	abort_message='Pb allocation variables couplees pour l''ecriture'
	1600	call abort_gcm(modname,abort_message,1)
	1601	endif
[138]	1602	endif
	1603
	1604	!
	1605	! Mise sur la bonne grille des champs a passer au coupleur
	1606	!
[147]	1607	cpl_index = 1
	1608	if (nisurf == is_sic) cpl_index = 2
	1609	call gath2cpl(cpl_sols(1,cpl_index), tmp_sols(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1610	call gath2cpl(cpl_nsol(1,cpl_index), tmp_nsol(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1611	call gath2cpl(cpl_rain(1,cpl_index), tmp_rain(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1612	call gath2cpl(cpl_snow(1,cpl_index), tmp_snow(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1613	call gath2cpl(cpl_evap(1,cpl_index), tmp_evap(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1614	call gath2cpl(cpl_tsol(1,cpl_index), tmp_tsol(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1615	call gath2cpl(cpl_fder(1,cpl_index), tmp_fder(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1616	call gath2cpl(cpl_albe(1,cpl_index), tmp_albe(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1617	call gath2cpl(cpl_taux(1,cpl_index), tmp_taux(1,1,cpl_index), klon, knon,iim,jjm, knindex)
	1618	call gath2cpl(cpl_tauy(1,cpl_index), tmp_tauy(1,1,cpl_index), klon, knon,iim,jjm, knindex)
[138]	1619
	1620	!
[105]	1621	! Si le domaine considere est la banquise, on envoie les champs au coupleur
	1622	!
[177]	1623	if (nisurf == is_sic .and. cumul) then
[105]	1624	wri_rain = 0.; wri_snow = 0.; wri_rcoa = 0.; wri_rriv = 0.
	1625	wri_taux = 0.; wri_tauy = 0.
	1626	call gath2cpl(pctsrf(1,is_oce), tamp_srf(1,1,1), klon, klon, iim, jjm, tamp_ind)
	1627	call gath2cpl(pctsrf(1,is_sic), tamp_srf(1,1,2), klon, klon, iim, jjm, tamp_ind)
[98]	1628
[105]	1629	wri_sol_ice = tmp_sols(:,:,2)
	1630	wri_sol_sea = tmp_sols(:,:,1)
	1631	wri_nsol_ice = tmp_nsol(:,:,2)
	1632	wri_nsol_sea = tmp_nsol(:,:,1)
	1633	wri_fder_ice = tmp_fder(:,:,2)
	1634	wri_evap_ice = tmp_evap(:,:,2)
	1635	wri_evap_sea = tmp_evap(:,:,1)
[290]	1636	!!$PB
	1637	wri_rriv = cpl_rriv(:,:)
	1638	wri_rcoa = cpl_rcoa(:,:)
[522]	1639	DO j = 1, jjm + 1
	1640	wri_calv(:,j) = sum(cpl_rlic(:,j)) / iim
	1641	enddo
[290]	1642
[105]	1643	where (tamp_zmasq /= 1.)
	1644	deno = tamp_srf(:,:,1) + tamp_srf(:,:,2)
	1645	wri_rain = tmp_rain(:,:,1) * tamp_srf(:,:,1) / deno + &
	1646	& tmp_rain(:,:,2) * tamp_srf(:,:,2) / deno
	1647	wri_snow = tmp_snow(:,:,1) * tamp_srf(:,:,1) / deno + &
	1648	& tmp_snow(:,:,2) * tamp_srf(:,:,2) / deno
	1649	wri_taux = tmp_taux(:,:,1) * tamp_srf(:,:,1) / deno + &
	1650	& tmp_taux(:,:,2) * tamp_srf(:,:,2) / deno
	1651	wri_tauy = tmp_tauy(:,:,1) * tamp_srf(:,:,1) / deno + &
	1652	& tmp_tauy(:,:,2) * tamp_srf(:,:,2) / deno
	1653	endwhere
[177]	1654	!
[313]	1655	! pour simuler la fonte des glaciers antarctiques
	1656	!
[394]	1657	!$$$ wri_rain = wri_rain &
	1658	!$$$ & + coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)
[522]	1659	! wri_calv = coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)
[313]	1660	!
[177]	1661	! on passe les coordonnées de la grille
	1662	!
[179]	1663
	1664	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,tmp_lon)
	1665	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,tmp_lat)
	1666
[177]	1667	DO i = 1, iim
	1668	tmp_lon(i,1) = rlon(i+1)
	1669	tmp_lon(i,jjm + 1) = rlon(i+1)
	1670	ENDDO
	1671	!
[179]	1672	! sortie netcdf des champs pour le changement de repere
	1673	!
	1674	ndexct(:)=0
[353]	1675	CALL histwrite(nidct,'tauxe',itau_w,wri_taux,iim*(jjm+1),ndexct)
	1676	CALL histwrite(nidct,'tauyn',itau_w,wri_tauy,iim*(jjm+1),ndexct)
	1677	CALL histwrite(nidct,'tmp_lon',itau_w,tmp_lon,iim*(jjm+1),ndexct)
	1678	CALL histwrite(nidct,'tmp_lat',itau_w,tmp_lat,iim*(jjm+1),ndexct)
[179]	1679
	1680	!
[177]	1681	! calcul 3 coordonnées du vent
	1682	!
	1683	CALL atm2geo (iim , jjm + 1, wri_taux, wri_tauy, tmp_lon, tmp_lat, &
	1684	& wri_tauxx, wri_tauyy, wri_tauzz )
[179]	1685	!
	1686	! sortie netcdf des champs apres changement de repere et juste avant
	1687	! envoi au coupleur
	1688	!
[353]	1689	CALL histwrite(nidct,cl_writ(1),itau_w,wri_sol_ice,iim*(jjm+1),ndexct)
	1690	CALL histwrite(nidct,cl_writ(2),itau_w,wri_sol_sea,iim*(jjm+1),ndexct)
	1691	CALL histwrite(nidct,cl_writ(3),itau_w,wri_nsol_ice,iim*(jjm+1),ndexct)
	1692	CALL histwrite(nidct,cl_writ(4),itau_w,wri_nsol_sea,iim*(jjm+1),ndexct)
	1693	CALL histwrite(nidct,cl_writ(5),itau_w,wri_fder_ice,iim*(jjm+1),ndexct)
	1694	CALL histwrite(nidct,cl_writ(6),itau_w,wri_evap_ice,iim*(jjm+1),ndexct)
	1695	CALL histwrite(nidct,cl_writ(7),itau_w,wri_evap_sea,iim*(jjm+1),ndexct)
	1696	CALL histwrite(nidct,cl_writ(8),itau_w,wri_rain,iim*(jjm+1),ndexct)
	1697	CALL histwrite(nidct,cl_writ(9),itau_w,wri_snow,iim*(jjm+1),ndexct)
	1698	CALL histwrite(nidct,cl_writ(10),itau_w,wri_rcoa,iim*(jjm+1),ndexct)
	1699	CALL histwrite(nidct,cl_writ(11),itau_w,wri_rriv,iim*(jjm+1),ndexct)
[394]	1700	CALL histwrite(nidct,cl_writ(12),itau_w,wri_calv,iim*(jjm+1),ndexct)
	1701	CALL histwrite(nidct,cl_writ(13),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
	1702	CALL histwrite(nidct,cl_writ(14),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
	1703	CALL histwrite(nidct,cl_writ(15),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
	1704	CALL histwrite(nidct,cl_writ(16),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
	1705	CALL histwrite(nidct,cl_writ(17),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
	1706	CALL histwrite(nidct,cl_writ(18),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
[179]	1707	CALL histsync(nidct)
	1708	! pas utile IF (lafin) CALL histclo(nidct)
[105]	1709	call intocpl(itime, (jjm+1)*iim, wri_sol_ice, wri_sol_sea, wri_nsol_ice,&
	1710	& wri_nsol_sea, wri_fder_ice, wri_evap_ice, wri_evap_sea, wri_rain, &
[394]	1711	& wri_snow, wri_rcoa, wri_rriv, wri_calv, wri_tauxx, wri_tauyy, &
	1712	& wri_tauzz, wri_tauxx, wri_tauyy, wri_tauzz,lafin )
[290]	1713	!
[105]	1714	cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
	1715	cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
[522]	1716	cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.; cpl_rlic = 0.
[105]	1717	!
	1718	! deallocation memoire variables temporaires
	1719	!
	1720	sum_error = 0
	1721	deallocate(tmp_sols, stat=error); sum_error = sum_error + error
	1722	deallocate(tmp_nsol, stat=error); sum_error = sum_error + error
	1723	deallocate(tmp_rain, stat=error); sum_error = sum_error + error
	1724	deallocate(tmp_snow, stat=error); sum_error = sum_error + error
	1725	deallocate(tmp_evap, stat=error); sum_error = sum_error + error
	1726	deallocate(tmp_fder, stat=error); sum_error = sum_error + error
	1727	deallocate(tmp_tsol, stat=error); sum_error = sum_error + error
	1728	deallocate(tmp_albe, stat=error); sum_error = sum_error + error
	1729	deallocate(tmp_taux, stat=error); sum_error = sum_error + error
	1730	deallocate(tmp_tauy, stat=error); sum_error = sum_error + error
[290]	1731	!!$PB
[398]	1732	!!$ deallocate(tmp_rriv, stat=error); sum_error = sum_error + error
	1733	!!$ deallocate(tmp_rcoa, stat=error); sum_error = sum_error + error
[105]	1734	if (sum_error /= 0) then
	1735	abort_message='Pb deallocation variables couplees'
	1736	call abort_gcm(modname,abort_message,1)
	1737	endif
	1738
	1739	endif
	1740
[138]	1741	endif ! fin (mod(itime, nexca) == 0)
[105]	1742	!
	1743	! on range les variables lues/sauvegardees dans les bonnes variables de sortie
	1744	!
	1745	if (nisurf == is_oce) then
[98]	1746	call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjm, knindex)
[105]	1747	else if (nisurf == is_sic) then
[140]	1748	call cpl2gath(read_sit, tsurf_new, klon, knon,iim,jjm, knindex)
	1749	call cpl2gath(read_alb_sic, alb_new, klon, knon,iim,jjm, knindex)
[98]	1750	endif
[140]	1751	pctsrf_new(:,nisurf) = pctsrf_sav(:,nisurf)
[115]	1752
	1753	! if (lafin) call quitcpl
[98]	1754
[90]	1755	END SUBROUTINE interfoce_cpl
	1756	!
	1757	!#########################################################################
	1758	!
[98]	1759
[90]	1760	SUBROUTINE interfoce_slab(nisurf)
	1761
	1762	! Cette routine sert d'interface entre le modele atmospherique et un
	1763	! modele de 'slab' ocean
	1764	!
	1765	! L. Fairhead 02/2000
	1766	!
	1767	! input:
	1768	! nisurf index de la surface a traiter (1 = sol continental)
	1769	!
	1770	! output:
	1771	!
	1772
	1773	! Parametres d'entree
	1774	integer, intent(IN) :: nisurf
	1775
	1776	END SUBROUTINE interfoce_slab
	1777	!
	1778	!#########################################################################
	1779	!
	1780	SUBROUTINE interfoce_lim(itime, dtime, jour, &
	1781	& klon, nisurf, knon, knindex, &
	1782	& debut, &
[109]	1783	& lmt_sst, pctsrf_new)
[90]	1784
	1785	! Cette routine sert d'interface entre le modele atmospherique et un fichier
	1786	! de conditions aux limites
	1787	!
	1788	! L. Fairhead 02/2000
	1789	!
	1790	! input:
	1791	! itime numero du pas de temps courant
	1792	! dtime pas de temps de la physique (en s)
	1793	! jour jour a lire dans l'annee
	1794	! nisurf index de la surface a traiter (1 = sol continental)
	1795	! knon nombre de points dans le domaine a traiter
	1796	! knindex index des points de la surface a traiter
	1797	! klon taille de la grille
	1798	! debut logical: 1er appel a la physique (initialisation)
	1799	!
	1800	! output:
	1801	! lmt_sst SST lues dans le fichier de CL
	1802	! pctsrf_new sous-maille fractionnelle
	1803	!
	1804
	1805
	1806	! Parametres d'entree
	1807	integer, intent(IN) :: itime
	1808	real , intent(IN) :: dtime
	1809	integer, intent(IN) :: jour
	1810	integer, intent(IN) :: nisurf
	1811	integer, intent(IN) :: knon
	1812	integer, intent(IN) :: klon
[147]	1813	integer, dimension(klon), intent(in) :: knindex
[90]	1814	logical, intent(IN) :: debut
	1815
	1816	! Parametres de sortie
[147]	1817	real, intent(out), dimension(klon) :: lmt_sst
[90]	1818	real, intent(out), dimension(klon,nbsrf) :: pctsrf_new
	1819
	1820	! Variables locales
	1821	integer :: ii
[112]	1822	INTEGER,save :: lmt_pas ! frequence de lecture des conditions limites
[90]	1823	! (en pas de physique)
	1824	logical,save :: deja_lu ! pour indiquer que le jour a lire a deja
	1825	! lu pour une surface precedente
	1826	integer,save :: jour_lu
	1827	integer :: ierr
	1828	character (len = 20) :: modname = 'interfoce_lim'
	1829	character (len = 80) :: abort_message
[179]	1830	character (len = 20),save :: fich ='limit.nc'
	1831	logical, save :: newlmt = .TRUE.
[295]	1832	logical, save :: check = .FALSE.
[90]	1833	! Champs lus dans le fichier de CL
[147]	1834	real, allocatable , save, dimension(:) :: sst_lu, rug_lu, nat_lu
[90]	1835	real, allocatable , save, dimension(:,:) :: pct_tmp
	1836	!
	1837	! quelques variables pour netcdf
	1838	!
	1839	#include "netcdf.inc"
	1840	integer :: nid, nvarid
	1841	integer, dimension(2) :: start, epais
	1842	!
	1843	! Fin déclaration
	1844	!
	1845
[112]	1846	if (debut .and. .not. allocated(sst_lu)) then
[90]	1847	lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
	1848	jour_lu = jour - 1
	1849	allocate(sst_lu(klon))
	1850	allocate(nat_lu(klon))
	1851	allocate(pct_tmp(klon,nbsrf))
	1852	endif
	1853
	1854	if ((jour - jour_lu) /= 0) deja_lu = .false.
	1855
[171]	1856	if (check) write(,)modname,' :: jour, jour_lu, deja_lu', jour, jour_lu, deja_lu
[140]	1857	if (check) write(,)modname,' :: itime, lmt_pas ', itime, lmt_pas,dtime
[90]	1858
	1859	! Tester d'abord si c'est le moment de lire le fichier
	1860	if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then
	1861	!
	1862	! Ouverture du fichier
	1863	!
[112]	1864	fich = trim(fich)
[90]	1865	ierr = NF_OPEN (fich, NF_NOWRITE,nid)
	1866	if (ierr.NE.NF_NOERR) then
	1867	abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
	1868	call abort_gcm(modname,abort_message,1)
	1869	endif
	1870	!
	1871	! La tranche de donnees a lire:
	1872	!
	1873	start(1) = 1
[478]	1874	start(2) = jour
[90]	1875	epais(1) = klon
	1876	epais(2) = 1
	1877	!
	1878	if (newlmt) then
	1879	!
	1880	! Fraction "ocean"
	1881	!
	1882	ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)
	1883	if (ierr /= NF_NOERR) then
	1884	abort_message = 'Le champ <FOCE> est absent'
	1885	call abort_gcm(modname,abort_message,1)
	1886	endif
	1887	#ifdef NC_DOUBLE
[112]	1888	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_oce))
[90]	1889	#else
[112]	1890	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_oce))
[90]	1891	#endif
	1892	if (ierr /= NF_NOERR) then
	1893	abort_message = 'Lecture echouee pour <FOCE>'
	1894	call abort_gcm(modname,abort_message,1)
	1895	endif
	1896	!
	1897	! Fraction "glace de mer"
	1898	!
	1899	ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)
	1900	if (ierr /= NF_NOERR) then
	1901	abort_message = 'Le champ <FSIC> est absent'
	1902	call abort_gcm(modname,abort_message,1)
	1903	endif
	1904	#ifdef NC_DOUBLE
[112]	1905	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_sic))
[90]	1906	#else
[112]	1907	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_sic))
[90]	1908	#endif
	1909	if (ierr /= NF_NOERR) then
	1910	abort_message = 'Lecture echouee pour <FSIC>'
	1911	call abort_gcm(modname,abort_message,1)
	1912	endif
	1913	!
	1914	! Fraction "terre"
	1915	!
	1916	ierr = NF_INQ_VARID(nid, 'FTER', nvarid)
	1917	if (ierr /= NF_NOERR) then
	1918	abort_message = 'Le champ <FTER> est absent'
	1919	call abort_gcm(modname,abort_message,1)
	1920	endif
	1921	#ifdef NC_DOUBLE
[112]	1922	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_ter))
[90]	1923	#else
[112]	1924	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_ter))
[90]	1925	#endif
	1926	if (ierr /= NF_NOERR) then
	1927	abort_message = 'Lecture echouee pour <FTER>'
	1928	call abort_gcm(modname,abort_message,1)
	1929	endif
	1930	!
	1931	! Fraction "glacier terre"
	1932	!
	1933	ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)
	1934	if (ierr /= NF_NOERR) then
	1935	abort_message = 'Le champ <FLIC> est absent'
	1936	call abort_gcm(modname,abort_message,1)
	1937	endif
	1938	#ifdef NC_DOUBLE
[112]	1939	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_lic))
[90]	1940	#else
[112]	1941	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_lic))
[90]	1942	#endif
	1943	if (ierr /= NF_NOERR) then
	1944	abort_message = 'Lecture echouee pour <FLIC>'
	1945	call abort_gcm(modname,abort_message,1)
	1946	endif
	1947	!
	1948	else ! on en est toujours a rnatur
	1949	!
	1950	ierr = NF_INQ_VARID(nid, 'NAT', nvarid)
	1951	if (ierr /= NF_NOERR) then
	1952	abort_message = 'Le champ <NAT> est absent'
	1953	call abort_gcm(modname,abort_message,1)
	1954	endif
	1955	#ifdef NC_DOUBLE
	1956	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, nat_lu)
	1957	#else
	1958	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, nat_lu)
	1959	#endif
	1960	if (ierr /= NF_NOERR) then
	1961	abort_message = 'Lecture echouee pour <NAT>'
	1962	call abort_gcm(modname,abort_message,1)
	1963	endif
	1964	!
	1965	! Remplissage des fractions de surface
	1966	! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice
	1967	!
	1968	pct_tmp = 0.0
	1969	do ii = 1, klon
	1970	pct_tmp(ii,nint(nat_lu(ii)) + 1) = 1.
	1971	enddo
	1972
	1973	!
	1974	! On se retrouve avec ocean en 1 et terre en 2 alors qu'on veut le contraire
	1975	!
	1976	pctsrf_new = pct_tmp
	1977	pctsrf_new (:,2)= pct_tmp (:,1)
	1978	pctsrf_new (:,1)= pct_tmp (:,2)
	1979	pct_tmp = pctsrf_new
	1980	endif ! fin test sur newlmt
	1981	!
	1982	! Lecture SST
	1983	!
	1984	ierr = NF_INQ_VARID(nid, 'SST', nvarid)
	1985	if (ierr /= NF_NOERR) then
	1986	abort_message = 'Le champ <SST> est absent'
	1987	call abort_gcm(modname,abort_message,1)
	1988	endif
	1989	#ifdef NC_DOUBLE
	1990	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, sst_lu)
	1991	#else
	1992	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, sst_lu)
	1993	#endif
	1994	if (ierr /= NF_NOERR) then
	1995	abort_message = 'Lecture echouee pour <SST>'
	1996	call abort_gcm(modname,abort_message,1)
	1997	endif
[109]	1998
[90]	1999	!
[109]	2000	! Fin de lecture
	2001	!
	2002	ierr = NF_CLOSE(nid)
	2003	deja_lu = .true.
	2004	jour_lu = jour
	2005	endif
	2006	!
	2007	! Recopie des variables dans les champs de sortie
	2008	!
[235]	2009	lmt_sst = 999999999.
[112]	2010	do ii = 1, knon
	2011	lmt_sst(ii) = sst_lu(knindex(ii))
	2012	enddo
[109]	2013
[171]	2014	pctsrf_new(:,is_oce) = pct_tmp(:,is_oce)
	2015	pctsrf_new(:,is_sic) = pct_tmp(:,is_sic)
	2016
[109]	2017	END SUBROUTINE interfoce_lim
	2018
	2019	!
	2020	!#########################################################################
	2021	!
	2022	SUBROUTINE interfsur_lim(itime, dtime, jour, &
	2023	& klon, nisurf, knon, knindex, &
	2024	& debut, &
	2025	& lmt_alb, lmt_rug)
	2026
	2027	! Cette routine sert d'interface entre le modele atmospherique et un fichier
	2028	! de conditions aux limites
	2029	!
	2030	! L. Fairhead 02/2000
	2031	!
	2032	! input:
	2033	! itime numero du pas de temps courant
	2034	! dtime pas de temps de la physique (en s)
	2035	! jour jour a lire dans l'annee
	2036	! nisurf index de la surface a traiter (1 = sol continental)
	2037	! knon nombre de points dans le domaine a traiter
	2038	! knindex index des points de la surface a traiter
	2039	! klon taille de la grille
	2040	! debut logical: 1er appel a la physique (initialisation)
	2041	!
	2042	! output:
	2043	! lmt_sst SST lues dans le fichier de CL
	2044	! lmt_alb Albedo lu
	2045	! lmt_rug longueur de rugosité lue
	2046	! pctsrf_new sous-maille fractionnelle
	2047	!
	2048
	2049
	2050	! Parametres d'entree
	2051	integer, intent(IN) :: itime
	2052	real , intent(IN) :: dtime
	2053	integer, intent(IN) :: jour
	2054	integer, intent(IN) :: nisurf
	2055	integer, intent(IN) :: knon
	2056	integer, intent(IN) :: klon
[147]	2057	integer, dimension(klon), intent(in) :: knindex
[109]	2058	logical, intent(IN) :: debut
	2059
	2060	! Parametres de sortie
[147]	2061	real, intent(out), dimension(klon) :: lmt_alb
	2062	real, intent(out), dimension(klon) :: lmt_rug
[109]	2063
	2064	! Variables locales
	2065	integer :: ii
[140]	2066	integer,save :: lmt_pas ! frequence de lecture des conditions limites
[109]	2067	! (en pas de physique)
	2068	logical,save :: deja_lu_sur! pour indiquer que le jour a lire a deja
	2069	! lu pour une surface precedente
	2070	integer,save :: jour_lu_sur
	2071	integer :: ierr
[121]	2072	character (len = 20) :: modname = 'interfsur_lim'
[109]	2073	character (len = 80) :: abort_message
[179]	2074	character (len = 20),save :: fich ='limit.nc'
	2075	logical,save :: newlmt = .false.
[458]	2076	logical,save :: check = .false.
[109]	2077	! Champs lus dans le fichier de CL
	2078	real, allocatable , save, dimension(:) :: alb_lu, rug_lu
	2079	!
	2080	! quelques variables pour netcdf
	2081	!
	2082	#include "netcdf.inc"
[179]	2083	integer ,save :: nid, nvarid
	2084	integer, dimension(2),save :: start, epais
[109]	2085	!
	2086	! Fin déclaration
	2087	!
	2088
	2089	if (debut) then
	2090	lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
	2091	jour_lu_sur = jour - 1
	2092	allocate(alb_lu(klon))
	2093	allocate(rug_lu(klon))
	2094	endif
	2095
[112]	2096	if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false.
[109]	2097
[112]	2098	if (check) write(,)modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, deja_lu_sur
[140]	2099	if (check) write(,)modname,':: itime, lmt_pas', itime, lmt_pas
[441]	2100	if (check) call flush(6)
[109]	2101
	2102	! Tester d'abord si c'est le moment de lire le fichier
	2103	if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu_sur) then
	2104	!
	2105	! Ouverture du fichier
	2106	!
[112]	2107	fich = trim(fich)
[441]	2108	IF (check) WRITE(,)modname,' ouverture fichier ',fich
	2109	if (check) CALL flush(6)
[109]	2110	ierr = NF_OPEN (fich, NF_NOWRITE,nid)
	2111	if (ierr.NE.NF_NOERR) then
	2112	abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
	2113	call abort_gcm(modname,abort_message,1)
	2114	endif
	2115	!
	2116	! La tranche de donnees a lire:
[112]	2117
[109]	2118	start(1) = 1
[479]	2119	start(2) = jour
[109]	2120	epais(1) = klon
	2121	epais(2) = 1
	2122	!
[90]	2123	! Lecture Albedo
	2124	!
	2125	ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
	2126	if (ierr /= NF_NOERR) then
	2127	abort_message = 'Le champ <ALB> est absent'
	2128	call abort_gcm(modname,abort_message,1)
	2129	endif
	2130	#ifdef NC_DOUBLE
	2131	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, alb_lu)
	2132	#else
	2133	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)
	2134	#endif
	2135	if (ierr /= NF_NOERR) then
	2136	abort_message = 'Lecture echouee pour <ALB>'
	2137	call abort_gcm(modname,abort_message,1)
	2138	endif
	2139	!
	2140	! Lecture rugosité
	2141	!
	2142	ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
	2143	if (ierr /= NF_NOERR) then
	2144	abort_message = 'Le champ <RUG> est absent'
	2145	call abort_gcm(modname,abort_message,1)
	2146	endif
	2147	#ifdef NC_DOUBLE
	2148	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, rug_lu)
	2149	#else
	2150	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)
	2151	#endif
	2152	if (ierr /= NF_NOERR) then
	2153	abort_message = 'Lecture echouee pour <RUG>'
	2154	call abort_gcm(modname,abort_message,1)
	2155	endif
	2156
	2157	!
	2158	! Fin de lecture
	2159	!
	2160	ierr = NF_CLOSE(nid)
[109]	2161	deja_lu_sur = .true.
	2162	jour_lu_sur = jour
[90]	2163	endif
	2164	!
	2165	! Recopie des variables dans les champs de sortie
	2166	!
[235]	2167	!!$ lmt_alb(:) = 0.0
	2168	!!$ lmt_rug(:) = 0.0
	2169	lmt_alb(:) = 999999.
	2170	lmt_rug(:) = 999999.
[112]	2171	DO ii = 1, knon
	2172	lmt_alb(ii) = alb_lu(knindex(ii))
	2173	lmt_rug(ii) = rug_lu(knindex(ii))
	2174	enddo
[90]	2175
[109]	2176	END SUBROUTINE interfsur_lim
[90]	2177
	2178	!
	2179	!#########################################################################
	2180	!
	2181
[147]	2182	SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime, &
[90]	2183	& tsurf, p1lay, cal, beta, coef1lay, ps, &
[441]	2184	& precip_rain, precip_snow, snow, qsurf, &
[90]	2185	& radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
	2186	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
[98]	2187	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
[90]	2188
	2189	! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
	2190	! une temperature de surface (au cas ou ok_veget = false)
	2191	!
	2192	! L. Fairhead 4/2000
	2193	!
	2194	! input:
	2195	! knon nombre de points a traiter
[98]	2196	! nisurf surface a traiter
[90]	2197	! tsurf temperature de surface
	2198	! p1lay pression 1er niveau (milieu de couche)
	2199	! cal capacite calorifique du sol
	2200	! beta evap reelle
	2201	! coef1lay coefficient d'echange
	2202	! ps pression au sol
[98]	2203	! precip_rain precipitations liquides
	2204	! precip_snow precipitations solides
	2205	! snow champs hauteur de neige
	2206	! runoff runoff en cas de trop plein
[90]	2207	! petAcoef coeff. A de la resolution de la CL pour t
	2208	! peqAcoef coeff. A de la resolution de la CL pour q
	2209	! petBcoef coeff. B de la resolution de la CL pour t
	2210	! peqBcoef coeff. B de la resolution de la CL pour q
	2211	! radsol rayonnement net aus sol (LW + SW)
	2212	! dif_grnd coeff. diffusion vers le sol profond
	2213	!
	2214	! output:
	2215	! tsurf_new temperature au sol
[441]	2216	! qsurf humidite de l'air au dessus du sol
[90]	2217	! fluxsens flux de chaleur sensible
	2218	! fluxlat flux de chaleur latente
	2219	! dflux_s derivee du flux de chaleur sensible / Ts
	2220	! dflux_l derivee du flux de chaleur latente / Ts
	2221	!
	2222
	2223	#include "YOETHF.inc"
	2224	#include "FCTTRE.inc"
[258]	2225	#include "indicesol.inc"
[90]	2226
	2227	! Parametres d'entree
[147]	2228	integer, intent(IN) :: knon, nisurf, klon
[90]	2229	real , intent(IN) :: dtime
[147]	2230	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
	2231	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
	2232	real, dimension(klon), intent(IN) :: ps, q1lay
	2233	real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
	2234	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
	2235	real, dimension(klon), intent(IN) :: radsol, dif_grnd
	2236	real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
[441]	2237	real, dimension(klon), intent(INOUT) :: snow, qsurf
[90]	2238
	2239	! Parametres sorties
[147]	2240	real, dimension(klon), intent(OUT):: tsurf_new, evap, fluxsens, fluxlat
	2241	real, dimension(klon), intent(OUT):: dflux_s, dflux_l
[90]	2242
	2243	! Variables locales
	2244	integer :: i
[147]	2245	real, dimension(klon) :: zx_mh, zx_nh, zx_oh
	2246	real, dimension(klon) :: zx_mq, zx_nq, zx_oq
	2247	real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
[157]	2248	real, dimension(klon) :: zx_sl, zx_k1
[181]	2249	real, dimension(klon) :: zx_q_0 , d_ts
[90]	2250	real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
[98]	2251	real :: bilan_f, fq_fonte
[177]	2252	REAL :: subli, fsno
[441]	2253	REAL :: qsat_new, q1_new
[90]	2254	real, parameter :: t_grnd = 271.35, t_coup = 273.15
[177]	2255	!! PB temporaire en attendant mieux pour le modele de neige
	2256	REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
	2257	!
[458]	2258	logical, save :: check = .false.
[90]	2259	character (len = 20) :: modname = 'calcul_fluxs'
[179]	2260	logical, save :: fonte_neige = .false.
	2261	real, save :: max_eau_sol = 150.0
[98]	2262	character (len = 80) :: abort_message
[179]	2263	logical,save :: first = .true.,second=.false.
[90]	2264
[140]	2265	if (check) write(,)'Entree ', modname,' surface = ',nisurf
[98]	2266
[441]	2267	IF (check) THEN
	2268	WRITE(,)' radsol (min, max)' &
	2269	& , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
	2270	CALL flush(6)
	2271	ENDIF
	2272
[275]	2273	if (size(coastalflow) /= knon .AND. nisurf == is_ter) then
[98]	2274	write(,)'Bizarre, le nombre de points continentaux'
	2275	write(,)'a change entre deux appels. J''arrete ...'
	2276	abort_message='Pb run_off'
	2277	call abort_gcm(modname,abort_message,1)
	2278	endif
	2279	!
	2280	! Traitement neige et humidite du sol
	2281	!
[258]	2282	!!$ WRITE(,)'test calcul_flux, surface ', nisurf
	2283	!!PB test
	2284	!!$ if (nisurf == is_oce) then
	2285	!!$ snow = 0.
	2286	!!$ qsol = max_eau_sol
	2287	!!$ else
	2288	!!$ where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
	2289	!!$ where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
	2290	!!$! snow = max(0.0, snow + (precip_snow - evap) * dtime)
	2291	!!$ where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
	2292	!!$ endif
[441]	2293	!!$ IF (nisurf /= is_ter) qsol = max_eau_sol
[98]	2294
	2295
[90]	2296	!
	2297	! Initialisation
	2298	!
[290]	2299	evap = 0.
	2300	fluxsens=0.
	2301	fluxlat=0.
	2302	dflux_s = 0.
	2303	dflux_l = 0.
[90]	2304	!
	2305	! zx_qs = qsat en kg/kg
	2306	!
	2307	DO i = 1, knon
	2308	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
	2309	IF (thermcep) THEN
	2310	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
	2311	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	2312	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
	2313	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
	2314	zx_qs=MIN(0.5,zx_qs)
	2315	zcor=1./(1.-retv*zx_qs)
	2316	zx_qs=zx_qs*zcor
	2317	zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
	2318	& /RLVTT / zx_pkh(i)
	2319	ELSE
	2320	IF (tsurf(i).LT.t_coup) THEN
	2321	zx_qs = qsats(tsurf(i)) / ps(i)
	2322	zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
	2323	& / zx_pkh(i)
	2324	ELSE
	2325	zx_qs = qsatl(tsurf(i)) / ps(i)
	2326	zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
	2327	& / zx_pkh(i)
	2328	ENDIF
	2329	ENDIF
	2330	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
	2331	zx_qsat(i) = zx_qs
	2332	zx_coef(i) = coef1lay(i) &
	2333	& * (1.0+SQRT(u1lay(i)2+v1lay(i)2)) &
	2334	& * p1lay(i)/(RD*t1lay(i))
[177]	2335
[90]	2336	ENDDO
	2337
	2338
	2339	! === Calcul de la temperature de surface ===
	2340	!
	2341	! zx_sl = chaleur latente d'evaporation ou de sublimation
	2342	!
	2343	do i = 1, knon
	2344	zx_sl(i) = RLVTT
	2345	if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
	2346	zx_k1(i) = zx_coef(i)
	2347	enddo
	2348
	2349
	2350	do i = 1, knon
	2351	! Q
	2352	zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
	2353	zx_mq(i) = beta(i) * zx_k1(i) * &
	2354	& (peqAcoef(i) - zx_qsat(i) &
	2355	& + zx_dq_s_dt(i) * tsurf(i)) &
	2356	& / zx_oq(i)
	2357	zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
	2358	& / zx_oq(i)
	2359
	2360	! H
	2361	zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
	2362	zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
	2363	zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
	2364
	2365	! Tsurface
	2366	tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
	2367	& (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
	2368	& + dif_grnd(i) * t_grnd * dtime)/ &
	2369	& ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
	2370	& zx_nh(i) + zx_sl(i) * zx_nq(i)) &
	2371	& + dtime * dif_grnd(i))
[98]	2372
	2373	!
	2374	! Y'a-t-il fonte de neige?
	2375	!
[181]	2376	! fonte_neige = (nisurf /= is_oce) .AND. &
	2377	! & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
	2378	! & .AND. (tsurf_new(i) >= RTT)
	2379	! if (fonte_neige) tsurf_new(i) = RTT
[90]	2380	d_ts(i) = tsurf_new(i) - tsurf(i)
[181]	2381	! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
	2382	! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
[90]	2383	!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
	2384	!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
[98]	2385	evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
	2386	fluxlat(i) = - evap(i) * zx_sl(i)
[90]	2387	fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
	2388	! Derives des flux dF/dTs (W m-2 K-1):
	2389	dflux_s(i) = zx_nh(i)
	2390	dflux_l(i) = (zx_sl(i) * zx_nq(i))
[441]	2391	! Nouvelle valeure de l'humidite au dessus du sol
	2392	qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
	2393	q1_new = peqAcoef(i) - peqBcoef(i)evap(i)dtime
	2394	qsurf(i)=q1_new(1.-beta(i)) + beta(i)qsat_new
[98]	2395	!
	2396	! en cas de fonte de neige
	2397	!
[181]	2398	! if (fonte_neige) then
	2399	! bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
	2400	! & dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
	2401	! & RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
	2402	! bilan_f = max(0., bilan_f)
	2403	! fq_fonte = bilan_f / zx_sl(i)
	2404	! snow(i) = max(0., snow(i) - fq_fonte * dtime)
	2405	! qsol(i) = qsol(i) + (fq_fonte * dtime)
	2406	! endif
[258]	2407	!!$ if (nisurf == is_ter) &
	2408	!!$ & run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
	2409	!!$ qsol(i) = min(qsol(i), max_eau_sol)
[90]	2410	ENDDO
	2411
	2412	END SUBROUTINE calcul_fluxs
	2413	!
	2414	!#########################################################################
	2415	!
[98]	2416	SUBROUTINE gath2cpl(champ_in, champ_out, klon, knon, iim, jjm, knindex)
	2417
	2418	! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
	2419	! au coupleur.
	2420	!
	2421	!
	2422	! input:
	2423	! champ_in champ sur la grille gathere
	2424	! knon nombre de points dans le domaine a traiter
	2425	! knindex index des points de la surface a traiter
	2426	! klon taille de la grille
	2427	! iim,jjm dimension de la grille 2D
	2428	!
	2429	! output:
	2430	! champ_out champ sur la grille 2D
	2431	!
	2432	! input
	2433	integer :: klon, knon, iim, jjm
[147]	2434	real, dimension(klon) :: champ_in
	2435	integer, dimension(klon) :: knindex
[98]	2436	! output
	2437	real, dimension(iim,jjm+1) :: champ_out
	2438	! local
	2439	integer :: i, ig, j
	2440	real, dimension(klon) :: tamp
	2441
[105]	2442	tamp = 0.
[98]	2443	do i = 1, knon
	2444	ig = knindex(i)
	2445	tamp(ig) = champ_in(i)
	2446	enddo
[140]	2447	ig = 1
	2448	champ_out(:,1) = tamp(ig)
[98]	2449	do j = 2, jjm
	2450	do i = 1, iim
[140]	2451	ig = ig + 1
	2452	champ_out(i,j) = tamp(ig)
[98]	2453	enddo
	2454	enddo
[140]	2455	ig = ig + 1
	2456	champ_out(:,jjm+1) = tamp(ig)
[98]	2457
	2458	END SUBROUTINE gath2cpl
	2459	!
	2460	!#########################################################################
	2461	!
	2462	SUBROUTINE cpl2gath(champ_in, champ_out, klon, knon, iim, jjm, knindex)
	2463
	2464	! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
	2465	! au coupleur.
	2466	!
	2467	!
	2468	! input:
	2469	! champ_in champ sur la grille gathere
	2470	! knon nombre de points dans le domaine a traiter
	2471	! knindex index des points de la surface a traiter
	2472	! klon taille de la grille
	2473	! iim,jjm dimension de la grille 2D
	2474	!
	2475	! output:
	2476	! champ_out champ sur la grille 2D
	2477	!
	2478	! input
	2479	integer :: klon, knon, iim, jjm
	2480	real, dimension(iim,jjm+1) :: champ_in
[147]	2481	integer, dimension(klon) :: knindex
[98]	2482	! output
[147]	2483	real, dimension(klon) :: champ_out
[98]	2484	! local
	2485	integer :: i, ig, j
	2486	real, dimension(klon) :: tamp
[179]	2487	logical ,save :: check = .false.
[98]	2488
[140]	2489	ig = 1
	2490	tamp(ig) = champ_in(1,1)
[98]	2491	do j = 2, jjm
	2492	do i = 1, iim
[140]	2493	ig = ig + 1
	2494	tamp(ig) = champ_in(i,j)
[98]	2495	enddo
	2496	enddo
[140]	2497	ig = ig + 1
	2498	tamp(ig) = champ_in(1,jjm+1)
[98]	2499
	2500	do i = 1, knon
	2501	ig = knindex(i)
	2502	champ_out(i) = tamp(ig)
	2503	enddo
	2504
	2505	END SUBROUTINE cpl2gath
	2506	!
	2507	!#########################################################################
	2508	!
[258]	2509	SUBROUTINE albsno(klon, knon,dtime,agesno,alb_neig_grid, precip_snow)
[109]	2510	IMPLICIT none
[112]	2511
[258]	2512	INTEGER :: klon, knon
[112]	2513	INTEGER, PARAMETER :: nvm = 8
[258]	2514	REAL :: dtime
[112]	2515	REAL, dimension(klon,nvm) :: veget
[258]	2516	REAL, DIMENSION(klon) :: alb_neig_grid, agesno, precip_snow
[112]	2517
	2518	INTEGER :: i, nv
	2519
	2520	REAL, DIMENSION(nvm),SAVE :: init, decay
	2521	REAL :: as
[109]	2522	DATA init /0.55, 0.14, 0.18, 0.29, 0.15, 0.15, 0.14, 0./
	2523	DATA decay/0.30, 0.67, 0.63, 0.45, 0.40, 0.14, 0.06, 1./
[112]	2524
[109]	2525	veget = 0.
	2526	veget(:,1) = 1. ! desert partout
[258]	2527	DO i = 1, knon
[112]	2528	alb_neig_grid(i) = 0.0
[109]	2529	ENDDO
	2530	DO nv = 1, nvm
[258]	2531	DO i = 1, knon
[109]	2532	as = init(nv)+decay(nv)*EXP(-agesno(i)/5.)
[112]	2533	alb_neig_grid(i) = alb_neig_grid(i) + veget(i,nv)*as
[109]	2534	ENDDO
	2535	ENDDO
[258]	2536	!
	2537	!! modilation en fonction de l'age de la neige
	2538	!
	2539	DO i = 1, knon
	2540	agesno(i) = (agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
	2541	& * EXP(-1.MAX(0.0,precip_snow(i))dtime/0.3)
	2542	agesno(i) = MAX(agesno(i),0.0)
	2543	ENDDO
[112]	2544
[109]	2545	END SUBROUTINE albsno
	2546	!
	2547	!#########################################################################
	2548	!
[181]	2549
	2550	SUBROUTINE fonte_neige( klon, knon, nisurf, dtime, &
	2551	& tsurf, p1lay, cal, beta, coef1lay, ps, &
	2552	& precip_rain, precip_snow, snow, qsol, &
	2553	& radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
	2554	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
[457]	2555	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
	2556	!IM cf JLD
[522]	2557	& fqcalving,ffonte,run_off_lic_0)
[181]	2558
	2559	! Routine de traitement de la fonte de la neige dans le cas du traitement
	2560	! de sol simplifié
	2561	!
	2562	! LF 03/2001
	2563	! input:
	2564	! knon nombre de points a traiter
	2565	! nisurf surface a traiter
	2566	! tsurf temperature de surface
	2567	! p1lay pression 1er niveau (milieu de couche)
	2568	! cal capacite calorifique du sol
	2569	! beta evap reelle
	2570	! coef1lay coefficient d'echange
	2571	! ps pression au sol
	2572	! precip_rain precipitations liquides
	2573	! precip_snow precipitations solides
	2574	! snow champs hauteur de neige
[441]	2575	! qsol hauteur d'eau contenu dans le sol
[181]	2576	! runoff runoff en cas de trop plein
	2577	! petAcoef coeff. A de la resolution de la CL pour t
	2578	! peqAcoef coeff. A de la resolution de la CL pour q
	2579	! petBcoef coeff. B de la resolution de la CL pour t
	2580	! peqBcoef coeff. B de la resolution de la CL pour q
	2581	! radsol rayonnement net aus sol (LW + SW)
	2582	! dif_grnd coeff. diffusion vers le sol profond
	2583	!
	2584	! output:
	2585	! tsurf_new temperature au sol
	2586	! fluxsens flux de chaleur sensible
	2587	! fluxlat flux de chaleur latente
	2588	! dflux_s derivee du flux de chaleur sensible / Ts
	2589	! dflux_l derivee du flux de chaleur latente / Ts
[522]	2590	! in/out:
	2591	! run_off_lic_0 run off glacier du pas de temps précedent
[181]	2592	!
	2593
	2594	#include "YOETHF.inc"
	2595	#include "FCTTRE.inc"
[258]	2596	#include "indicesol.inc"
[457]	2597	!IM cf JLD
[522]	2598	#include "YOMCST.inc"
[181]	2599
	2600	! Parametres d'entree
	2601	integer, intent(IN) :: knon, nisurf, klon
	2602	real , intent(IN) :: dtime
	2603	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
	2604	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
	2605	real, dimension(klon), intent(IN) :: ps, q1lay
	2606	real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
	2607	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
	2608	real, dimension(klon), intent(IN) :: radsol, dif_grnd
	2609	real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
	2610	real, dimension(klon), intent(INOUT) :: snow, qsol
	2611
	2612	! Parametres sorties
	2613	real, dimension(klon), intent(INOUT):: tsurf_new, evap, fluxsens, fluxlat
	2614	real, dimension(klon), intent(INOUT):: dflux_s, dflux_l
[457]	2615	! Flux thermique utiliser pour fondre la neige
	2616	real, dimension(klon), intent(INOUT):: ffonte
	2617	! Flux d'eau "perdue" par la surface et necessaire pour que limiter la
	2618	! hauteur de neige, en kg/m2/s
	2619	real, dimension(klon), intent(INOUT):: fqcalving
[522]	2620	real, dimension(klon), intent(INOUT):: run_off_lic_0
[181]	2621	! Variables locales
[457]	2622	! Masse maximum de neige (kg/m2). Au dessus de ce seuil, la neige
	2623	! en exces "s'ecoule" (calving)
[486]	2624	! real, parameter :: snow_max=1.
	2625	!IM cf JLD/GK
	2626	real, parameter :: snow_max=3000.
[181]	2627	integer :: i
	2628	real, dimension(klon) :: zx_mh, zx_nh, zx_oh
	2629	real, dimension(klon) :: zx_mq, zx_nq, zx_oq
	2630	real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
	2631	real, dimension(klon) :: zx_sl, zx_k1
	2632	real, dimension(klon) :: zx_q_0 , d_ts
	2633	real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
	2634	real :: bilan_f, fq_fonte
	2635	REAL :: subli, fsno
[522]	2636	REAL, DIMENSION(klon) :: bil_eau_s, snow_evap
[181]	2637	real, parameter :: t_grnd = 271.35, t_coup = 273.15
	2638	!! PB temporaire en attendant mieux pour le modele de neige
[457]	2639	! REAL, parameter :: chasno = RLMLT/(2.3867E+06*0.15)
[181]	2640	REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
[486]	2641	!IM cf JLD/ GKtest
	2642	REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)
	2643	! fin GKtest
[181]	2644	!
[295]	2645	logical, save :: check = .FALSE.
[181]	2646	character (len = 20) :: modname = 'fonte_neige'
	2647	logical, save :: neige_fond = .false.
	2648	real, save :: max_eau_sol = 150.0
	2649	character (len = 80) :: abort_message
	2650	logical,save :: first = .true.,second=.false.
[522]	2651	real :: coeff_rel
[181]	2652
[522]	2653
[181]	2654	if (check) write(,)'Entree ', modname,' surface = ',nisurf
	2655
	2656	! Initialisations
[522]	2657	coeff_rel = dtime/(tau_calv * rday)
[441]	2658	bil_eau_s(:) = 0.
[181]	2659	DO i = 1, knon
	2660	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
	2661	IF (thermcep) THEN
	2662	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
	2663	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	2664	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
	2665	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
	2666	zx_qs=MIN(0.5,zx_qs)
	2667	zcor=1./(1.-retv*zx_qs)
	2668	zx_qs=zx_qs*zcor
	2669	zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
	2670	& /RLVTT / zx_pkh(i)
	2671	ELSE
	2672	IF (tsurf(i).LT.t_coup) THEN
	2673	zx_qs = qsats(tsurf(i)) / ps(i)
	2674	zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
	2675	& / zx_pkh(i)
	2676	ELSE
	2677	zx_qs = qsatl(tsurf(i)) / ps(i)
	2678	zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
	2679	& / zx_pkh(i)
	2680	ENDIF
	2681	ENDIF
	2682	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
	2683	zx_qsat(i) = zx_qs
	2684	zx_coef(i) = coef1lay(i) &
	2685	& * (1.0+SQRT(u1lay(i)2+v1lay(i)2)) &
	2686	& * p1lay(i)/(RD*t1lay(i))
	2687	ENDDO
[223]	2688
	2689
[181]	2690	! === Calcul de la temperature de surface ===
	2691	!
	2692	! zx_sl = chaleur latente d'evaporation ou de sublimation
	2693	!
	2694	do i = 1, knon
	2695	zx_sl(i) = RLVTT
	2696	if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
	2697	zx_k1(i) = zx_coef(i)
	2698	enddo
	2699
	2700
	2701	do i = 1, knon
	2702	! Q
	2703	zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
	2704	zx_mq(i) = beta(i) * zx_k1(i) * &
	2705	& (peqAcoef(i) - zx_qsat(i) &
	2706	& + zx_dq_s_dt(i) * tsurf(i)) &
	2707	& / zx_oq(i)
	2708	zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
	2709	& / zx_oq(i)
	2710
	2711	! H
	2712	zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
	2713	zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
	2714	zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
	2715	enddo
	2716
[258]	2717
	2718	WHERE (precip_snow > 0.) snow = snow + (precip_snow * dtime)
[522]	2719	snow_evap = 0.
	2720	WHERE (evap > 0. )
	2721	snow_evap = MIN (snow / dtime, evap)
	2722	snow = snow - snow_evap * dtime
	2723	snow = MAX(0.0, snow)
	2724	end where
	2725
	2726	! bil_eau_s = bil_eau_s + (precip_rain * dtime) - (evap - snow_evap) * dtime
	2727	bil_eau_s = (precip_rain * dtime) - (evap - snow_evap) * dtime
	2728
[181]	2729	!
	2730	! Y'a-t-il fonte de neige?
	2731	!
[457]	2732	ffonte=0.
[181]	2733	do i = 1, knon
	2734	neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
	2735	& .AND. tsurf_new(i) >= RTT)
	2736	if (neige_fond) then
[258]	2737	fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i))
[457]	2738	ffonte(i) = fq_fonte * RLMLT/dtime
[258]	2739	snow(i) = max(0., snow(i) - fq_fonte)
[441]	2740	bil_eau_s(i) = bil_eau_s(i) + fq_fonte
[258]	2741	tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno
[457]	2742	!IM cf JLD OK
[486]	2743	!IM cf JLD/ GKtest fonte aussi pour la glace
	2744	IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
	2745	fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,0.0)
	2746	ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime
	2747	bil_eau_s(i) = bil_eau_s(i) + fq_fonte
	2748	tsurf_new(i) = RTT
	2749	ENDIF
[181]	2750	d_ts(i) = tsurf_new(i) - tsurf(i)
	2751	endif
[457]	2752	!
	2753	! s'il y a une hauteur trop importante de neige, elle s'coule
	2754	fqcalving(i) = max(0., snow(i) - snow_max)/dtime
	2755	snow(i)=min(snow(i),snow_max)
	2756	!
[441]	2757	IF (nisurf == is_ter) then
	2758	qsol(i) = qsol(i) + bil_eau_s(i)
	2759	run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)
	2760	qsol(i) = MIN(qsol(i), max_eau_sol)
[522]	2761	else if (nisurf == is_lic) then
	2762	run_off_lic(i) = (coeff_rel * fqcalving(i)) + &
	2763	& (1. - coeff_rel) * run_off_lic_0(i)
	2764	run_off_lic_0(i) = run_off_lic(i)
	2765	run_off_lic(i) = run_off_lic(i) + bil_eau_s(i)/dtime
[441]	2766	endif
[181]	2767	enddo
	2768
	2769	END SUBROUTINE fonte_neige
	2770	!
	2771	!#########################################################################
	2772	!
[90]	2773	END MODULE interface_surf

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

Download in other formats:

Original Format