Context Navigation

interface_surf.F90 @ 544

Last change on this file since 544 was 544, checked in by lmdzadmin, 20 years ago
Incorporation des modifications necessaires a l'utilisation de la librairie Psmile/PRISM, et creation d'un tag IPSL-CM4_PSMILE, selon M.-E. Demory LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 95.5 KB

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

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

Download in other formats:

Original Format