Context Navigation

interface_surf.F90 @ 365

Last change on this file since 365 was 365, checked in by lmdzadmin, 22 years ago
L'humidite saturante de surface calculee par ORCHIDEE n'etait pas repassee a clmain. IM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 88.7 KB

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

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

Download in other formats:

Original Format