Context Navigation

interface_surf.F90 @ 522

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

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

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

Download in other formats:

Original Format