Context Navigation

interface_surf.F90 @ 557

Last change on this file since 557 was 557, checked in by lmdzadmin, 20 years ago
Initialisations diverses heritees de LMDZ.3.3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 95.1 KB

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

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

Download in other formats:

Original Format