Context Navigation

interface_surf.F90 @ 224

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

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

Download in other formats:

Original Format