Context Navigation

interface_surf.F90 @ 398

Last change on this file since 398 was 398, checked in by lmdzadmin, 22 years ago
Réglage decalage runoff entre orchidee et coupleur OM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 89.0 KB

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

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

Download in other formats:

Original Format