Context Navigation

interface_surf.F90 @ 295

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

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

Download in other formats:

Original Format