Context Navigation

interface_surf.F90 @ 731

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

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

Download in other formats:

Original Format