Context Navigation

interface_surf.F90 @ 594

Last change on this file since 594 was 594, checked in by Laurent Fairhead, 19 years ago
Bugs introduit par la synchro avec LOOP je sais pas comment GG LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 97.3 KB

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

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

Download in other formats:

Original Format