Context Navigation

interface_surf.F90 @ 546

Last change on this file since 546 was 546, checked in by lmdzadmin, 20 years ago
Suite des adaptations au portage sur SGI et VPP pour Prism (Caubel & Demory) LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 95.5 KB

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

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

Download in other formats:

Original Format