Context Navigation

interface_surf.F90 @ 461

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

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

Download in other formats:

Original Format