Context Navigation

interface_surf.F90 @ 665

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

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

Download in other formats:

Original Format