Context Navigation

interface_surf.F90 @ 5000

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

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

Download in other formats:

Original Format