Context Navigation

interface_surf.F90 @ 467

Last change on this file since 467 was 467, checked in by lmdzadmin, 21 years ago
Modifs sur les seuils (cdrag etc...), inclusion des diagnostics ISCCP par Ionela LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 94.0 KB

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

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

Download in other formats:

Original Format