Context Navigation

interface_surf.F90 @ 232

Last change on this file since 232 was 229, checked in by lmdzadmin, 23 years ago
Merge avec branche principale LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 95.7 KB

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

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

Download in other formats:

Original Format