Context Navigation

interface_surf.F90 @ 441

Last change on this file since 441 was 441, checked in by lmdzadmin, 21 years ago
Remplace qsol par qsurf JLD Rajout zlev, swdown_vrai, albedo_keep JP IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 90.5 KB

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

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

Download in other formats:

Original Format