Context Navigation

interface_surf.F90 @ 1097

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

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

Download in other formats:

Original Format