Context Navigation

interface_surf.F90 @ 314

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

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

Download in other formats:

Original Format