Context Navigation

interface_surf.F90 @ 1308

Last change on this file since 1308 was 564, checked in by lmdzadmin, 20 years ago
Encore quelques initialisations pour le portage AC LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 95.5 KB

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

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

Download in other formats:

Original Format