Context Navigation

interface_surf.F90 @ 520

Last change on this file since 520 was 513, checked in by lmdzadmin, 21 years ago
Initialisation albedo lic a 0.77 IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 94.8 KB

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

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

Download in other formats:

Original Format