Context Navigation

interface_surf.F90 @ 258

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

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

Download in other formats:

Original Format