Context Navigation

interface_surf.F90 @ 353

Last change on this file since 353 was 353, checked in by lmdzadmin, 23 years ago

2 changements pour les fichiers histoire:

utilisation de l'entree "rectilineaire" de IOIPSL pour ne plus avoir

lancer ncregular a chaque fois

le calendrier des fichiers histoire est maintenant base sur la date d'initialisation de la simulation plutot que sur la date de depart du

job

en cours

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format