Context Navigation

interface_surf.F90 @ 590

Last change on this file since 590 was 590, checked in by Laurent Fairhead, 19 years ago
Synchro avec LOOP PC LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 97.3 KB

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

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

Download in other formats:

Original Format