Context Navigation

interface_surf.F90 @ 3414

Last change on this file since 3414 was 716, checked in by Laurent Fairhead, 18 years ago
Petits probleme de declaration LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 95.5 KB

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

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

Download in other formats:

Original Format