Context Navigation

interface_surf.F90 @ 589

Last change on this file since 589 was 560, checked in by lmdzadmin, 20 years ago
pb de commentaires LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 95.1 KB

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

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

Download in other formats:

Original Format