Context Navigation

interface_surf.F90 @ 623

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

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

Download in other formats:

Original Format