Context Navigation

interface_surf.F90 @ 241

Last change on this file since 241 was 236, checked in by lmdzadmin, 24 years ago
Inclus nouveau soil dans interface_surf LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 85.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	zfra = MAX(0.0,MIN(1.0,snow/(snow+10.0)))
348	alb_new = alb_neigzfra + alb_new(1.0-zfra)
349	z0_new = SQRT(z0_new2+rugoro2)
350
351	else
352	CALL albsno(klon,agesno,alb_neig_grid)
353	!
354	! appel a sechiba
355	!
356	call interfsol(itime, klon, dtime, date0, nisurf, knon, &
357	& knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
358	& debut, lafin, ok_veget, &
359	& zlev, u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
360	& tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
361	& precip_rain, precip_snow, sollwdown, swnet, swdown, &
362	& tsurf, p1lay/100., ps, radsol, &
363	& evap, fluxsens, fluxlat, &
364	& tsol_rad, tsurf_new, alb_new, emis_new, z0_new, dflux_l, dflux_s)
365
366	!
367	! ajout de la contribution du relief
368	!
369	z0_new = SQRT(z0_new2+rugoro2)
370
371	endif
372	!
373	! Remplissage des pourcentages de surface
374	!
375	pctsrf_new(:,nisurf) = pctsrf(:,nisurf)
376
377	else if (nisurf == is_oce) then
378
379	if (check) write(,)'ocean, nisurf = ',nisurf
380
381
382	!
383	! Surface "ocean" appel a l'interface avec l'ocean
384	!
385	if (ocean == 'couple') then
386	! nexca = 0
387	if (nexca == 0) then
388	abort_message='nexca = 0 dans interfoce_cpl'
389	call abort_gcm(modname,abort_message,1)
390	endif
391
392	cumul = .false.
393
394	call interfoce(itime, dtime, cumul, &
395	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
396	& ocean, npas, nexca, debut, lafin, &
397	& swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
398	& fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
399	& tsurf_new, alb_new, alb_ice, pctsrf_new)
400
401	! else if (ocean == 'slab ') then
402	! call interfoce(nisurf)
403	else ! lecture conditions limites
404	call interfoce(itime, dtime, jour, &
405	& klon, nisurf, knon, knindex, &
406	& debut, &
407	& tsurf_new, pctsrf_new)
408
409	endif
410
411	tsurf_temp = tsurf_new
412	cal = 0.
413	beta = 1.
414	dif_grnd = 0.
415
416	call calcul_fluxs( klon, knon, nisurf, dtime, &
417	& tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
418	& precip_rain, precip_snow, snow, qsol, &
419	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
420	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
421	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
422
423	fder = fder + dflux_s + dflux_l
424
425	!
426	! 2eme appel a interfoce pour le cumul des champs (en particulier
427	! fluxsens et fluxlat calcules dans calcul_fluxs)
428	!
429	if (ocean == 'couple') then
430
431	cumul = .true.
432
433	call interfoce(itime, dtime, cumul, &
434	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
435	& ocean, npas, nexca, debut, lafin, &
436	& swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
437	& fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
438	& tsurf_new, alb_new, alb_ice, pctsrf_new)
439
440	! else if (ocean == 'slab ') then
441	! call interfoce(nisurf)
442
443	endif
444
445	!
446	! calcul albedo
447	!
448
449	if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then
450	CALL alboc(FLOAT(jour),rlat,alb_eau)
451	else ! cycle diurne
452	CALL alboc_cd(rmu0,alb_eau)
453	endif
454	DO ii =1, knon
455	alb_new(ii) = alb_eau(knindex(ii))
456	enddo
457
458	z0_new = rugos
459	!
460	else if (nisurf == is_sic) then
461
462	if (check) write(,)'sea ice, nisurf = ',nisurf
463
464	!
465	! Surface "glace de mer" appel a l'interface avec l'ocean
466	!
467	!
468	if (ocean == 'couple') then
469
470	cumul =.false.
471
472	call interfoce(itime, dtime, cumul, &
473	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
474	& ocean, npas, nexca, debut, lafin, &
475	& swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
476	& fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
477	& tsurf_new, alb_new, alb_ice, pctsrf_new)
478
479	tsurf_temp = tsurf_new
480	cal = 0.
481	dif_grnd = 0.
482	beta = 1.0
483
484	! else if (ocean == 'slab ') then
485	! call interfoce(nisurf)
486	ELSE
487	! ! lecture conditions limites
488	CALL interfoce(itime, dtime, jour, &
489	& klon, nisurf, knon, knindex, &
490	& debut, &
491	& tsurf_new, pctsrf_new)
492
493	CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
494
495	IF (soil_model) THEN
496	CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)
497	cal(1:knon) = RCPD / soilcap(1:knon)
498	radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
499	dif_grnd = 0.
500	!!$ WRITE(,) 'radsol'
501	!!$ WRITE(,) radsol(1 : knon)
502	!!$ WRITE(,) 'soilflux'
503	!!$ WRITE(,) soilflux(1 : knon)
504	ELSE
505	! if (check) write(,)'Sortie calbeta'
506	! if (check) write(,)'RCPD = ',RCPD,' capsol = '
507	! if (check) write(,)capsol
508	dif_grnd = 1.0 / tau_gl
509	cal = RCPD * calice
510	WHERE (snow > 0.0) cal = RCPD * calsno
511	ENDIF
512	tsurf_temp = tsurf
513	beta = 1.0
514	ENDIF
515
516	CALL calcul_fluxs( klon, knon, nisurf, dtime, &
517	& tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
518	& precip_rain, precip_snow, snow, qsol, &
519	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
520	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
521	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
522
523	IF (ocean /= 'couple') THEN
524	CALL fonte_neige( klon, knon, nisurf, dtime, &
525	& tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
526	& precip_rain, precip_snow, snow, qsol, &
527	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
528	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
529	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
530	ENDIF
531	!
532	! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
533	!
534	if (ocean == 'couple') then
535
536	cumul =.true.
537
538	call interfoce(itime, dtime, cumul, &
539	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
540	& ocean, npas, nexca, debut, lafin, &
541	& swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
542	& fluxlat, fluxsens, fder, albedo, taux, tauy, zmasq, &
543	& tsurf_new, alb_new, alb_ice, pctsrf_new)
544
545	! else if (ocean == 'slab ') then
546	! call interfoce(nisurf)
547
548	endif
549
550	!
551	! calcul albedo
552	!
553	!
554	zfra = MAX(0.0,MIN(1.0,snow/(snow+10.0)))
555	DO ii = 1, knon
556	index = knindex(ii)
557	alb_neig(ii) = alb_neig_grid(index)
558	ENDDO
559	alb_new = alb_neigzfra + 0.6 (1.0-zfra)
560
561	z0_new = 0.001
562
563	else if (nisurf == is_lic) then
564
565	if (check) write(,)'glacier, nisurf = ',nisurf
566
567	!
568	! Surface "glacier continentaux" appel a l'interface avec le sol
569	!
570	! call interfsol(nisurf)
571	IF (soil_model) THEN
572	CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil,soilcap, soilflux)
573	cal(1:knon) = RCPD / soilcap(1:knon)
574	radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
575	!!$ WRITE(,) 'radsol'
576	!!$ WRITE(*,'(16f17.4)') radsol(1 : knon)
577	!!$ WRITE(,) 'soilflux'
578	!!$ WRITE(*,'(16f17.4)')soilflux(1:knon)
579	ELSE
580	cal = RCPD * calice
581	WHERE (snow > 0.0) cal = RCPD * calsno
582	ENDIF
583	beta = 1.0
584	dif_grnd = 0.0
585
586	call calcul_fluxs( klon, knon, nisurf, dtime, &
587	& tsurf, p1lay, cal, beta, tq_cdrag, ps, &
588	& precip_rain, precip_snow, snow, qsol, &
589	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
590	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
591	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
592
593	call fonte_neige( klon, knon, nisurf, dtime, &
594	& tsurf, p1lay, cal, beta, tq_cdrag, ps, &
595	& precip_rain, precip_snow, snow, qsol, &
596	& radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
597	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
598	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
599
600	!
601	! calcul albedo
602	!
603	zfra = MAX(0.0,MIN(1.0,snow/(snow+10.0)))
604	DO ii = 1, knon
605	index = knindex(ii)
606	alb_neig(ii) = alb_neig_grid(index)
607	agesno(index) = (agesno(index) + (1.-agesno(index)/50.)*dtime/86400.)&
608	& * EXP(-1.MAX(0.0,precip_snow(ii))dtime/0.3)
609	agesno(index) = MAX(agesno(index),0.0)
610	IF(snow(ii) .LT. 0.0001) agesno(index) = 0.
611	ENDDO
612	alb_new = alb_neigzfra + 0.6 (1.0-zfra)
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
2000	! Parametres d'entree
2001	integer, intent(IN) :: knon, nisurf, klon
2002	real , intent(IN) :: dtime
2003	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
2004	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
2005	real, dimension(klon), intent(IN) :: ps, q1lay
2006	real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
2007	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
2008	real, dimension(klon), intent(IN) :: radsol, dif_grnd
2009	real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
2010	real, dimension(klon), intent(INOUT) :: snow, qsol
2011
2012	! Parametres sorties
2013	real, dimension(klon), intent(OUT):: tsurf_new, evap, fluxsens, fluxlat
2014	real, dimension(klon), intent(OUT):: dflux_s, dflux_l
2015
2016	! Variables locales
2017	integer :: i
2018	real, dimension(klon) :: zx_mh, zx_nh, zx_oh
2019	real, dimension(klon) :: zx_mq, zx_nq, zx_oq
2020	real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
2021	real, dimension(klon) :: zx_sl, zx_k1
2022	real, dimension(klon) :: zx_q_0 , d_ts
2023	real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
2024	real :: bilan_f, fq_fonte
2025	REAL :: subli, fsno
2026	real, parameter :: t_grnd = 271.35, t_coup = 273.15
2027	!! PB temporaire en attendant mieux pour le modele de neige
2028	REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
2029	!
2030	logical, save :: check = .true.
2031	character (len = 20) :: modname = 'calcul_fluxs'
2032	logical, save :: fonte_neige = .false.
2033	real, save :: max_eau_sol = 150.0
2034	character (len = 80) :: abort_message
2035	logical,save :: first = .true.,second=.false.
2036
2037	if (check) write(,)'Entree ', modname,' surface = ',nisurf
2038
2039	if (size(run_off) /= knon .AND. nisurf == is_ter) then
2040	write(,)'Bizarre, le nombre de points continentaux'
2041	write(,)'a change entre deux appels. J''arrete ...'
2042	abort_message='Pb run_off'
2043	call abort_gcm(modname,abort_message,1)
2044	endif
2045	!
2046	! Traitement neige et humidite du sol
2047	!
2048	if (nisurf == is_oce) then
2049	snow = 0.
2050	qsol = max_eau_sol
2051	else
2052	snow = snow + (precip_snow * dtime)
2053	where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
2054	! snow = max(0.0, snow + (precip_snow - evap) * dtime)
2055	qsol = qsol + (precip_rain - evap) * dtime
2056	endif
2057	IF (nisurf /= is_ter) qsol = max_eau_sol
2058
2059
2060	!
2061	! Initialisation
2062	!
2063	!
2064	! zx_qs = qsat en kg/kg
2065	!
2066	DO i = 1, knon
2067	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
2068	IF (thermcep) THEN
2069	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
2070	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
2071	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
2072	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
2073	zx_qs=MIN(0.5,zx_qs)
2074	zcor=1./(1.-retv*zx_qs)
2075	zx_qs=zx_qs*zcor
2076	zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
2077	& /RLVTT / zx_pkh(i)
2078	ELSE
2079	IF (tsurf(i).LT.t_coup) THEN
2080	zx_qs = qsats(tsurf(i)) / ps(i)
2081	zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
2082	& / zx_pkh(i)
2083	ELSE
2084	zx_qs = qsatl(tsurf(i)) / ps(i)
2085	zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
2086	& / zx_pkh(i)
2087	ENDIF
2088	ENDIF
2089	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
2090	zx_qsat(i) = zx_qs
2091	zx_coef(i) = coef1lay(i) &
2092	& * (1.0+SQRT(u1lay(i)2+v1lay(i)2)) &
2093	& * p1lay(i)/(RD*t1lay(i))
2094
2095	ENDDO
2096
2097
2098	! === Calcul de la temperature de surface ===
2099	!
2100	! zx_sl = chaleur latente d'evaporation ou de sublimation
2101	!
2102	do i = 1, knon
2103	zx_sl(i) = RLVTT
2104	if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
2105	zx_k1(i) = zx_coef(i)
2106	enddo
2107
2108
2109	do i = 1, knon
2110	! Q
2111	zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
2112	zx_mq(i) = beta(i) * zx_k1(i) * &
2113	& (peqAcoef(i) - zx_qsat(i) &
2114	& + zx_dq_s_dt(i) * tsurf(i)) &
2115	& / zx_oq(i)
2116	zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
2117	& / zx_oq(i)
2118
2119	! H
2120	zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
2121	zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
2122	zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
2123
2124	! Tsurface
2125	tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
2126	& (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
2127	& + dif_grnd(i) * t_grnd * dtime)/ &
2128	& ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
2129	& zx_nh(i) + zx_sl(i) * zx_nq(i)) &
2130	& + dtime * dif_grnd(i))
2131
2132	!
2133	! Y'a-t-il fonte de neige?
2134	!
2135	! fonte_neige = (nisurf /= is_oce) .AND. &
2136	! & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
2137	! & .AND. (tsurf_new(i) >= RTT)
2138	! if (fonte_neige) tsurf_new(i) = RTT
2139	d_ts(i) = tsurf_new(i) - tsurf(i)
2140	! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
2141	! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
2142	!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
2143	!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
2144	evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
2145	fluxlat(i) = - evap(i) * zx_sl(i)
2146	fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
2147	! Derives des flux dF/dTs (W m-2 K-1):
2148	dflux_s(i) = zx_nh(i)
2149	dflux_l(i) = (zx_sl(i) * zx_nq(i))
2150
2151	!
2152	! en cas de fonte de neige
2153	!
2154	! if (fonte_neige) then
2155	! bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
2156	! & dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
2157	! & RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
2158	! bilan_f = max(0., bilan_f)
2159	! fq_fonte = bilan_f / zx_sl(i)
2160	! snow(i) = max(0., snow(i) - fq_fonte * dtime)
2161	! qsol(i) = qsol(i) + (fq_fonte * dtime)
2162	! endif
2163	if (nisurf == is_ter) &
2164	& run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
2165	qsol(i) = min(qsol(i), max_eau_sol)
2166	ENDDO
2167
2168	END SUBROUTINE calcul_fluxs
2169	!
2170	!#########################################################################
2171	!
2172	SUBROUTINE gath2cpl(champ_in, champ_out, klon, knon, iim, jjm, knindex)
2173
2174	! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
2175	! au coupleur.
2176	!
2177	!
2178	! input:
2179	! champ_in champ sur la grille gathere
2180	! knon nombre de points dans le domaine a traiter
2181	! knindex index des points de la surface a traiter
2182	! klon taille de la grille
2183	! iim,jjm dimension de la grille 2D
2184	!
2185	! output:
2186	! champ_out champ sur la grille 2D
2187	!
2188	! input
2189	integer :: klon, knon, iim, jjm
2190	real, dimension(klon) :: champ_in
2191	integer, dimension(klon) :: knindex
2192	! output
2193	real, dimension(iim,jjm+1) :: champ_out
2194	! local
2195	integer :: i, ig, j
2196	real, dimension(klon) :: tamp
2197
2198	tamp = 0.
2199	do i = 1, knon
2200	ig = knindex(i)
2201	tamp(ig) = champ_in(i)
2202	enddo
2203	ig = 1
2204	champ_out(:,1) = tamp(ig)
2205	do j = 2, jjm
2206	do i = 1, iim
2207	ig = ig + 1
2208	champ_out(i,j) = tamp(ig)
2209	enddo
2210	enddo
2211	ig = ig + 1
2212	champ_out(:,jjm+1) = tamp(ig)
2213
2214	END SUBROUTINE gath2cpl
2215	!
2216	!#########################################################################
2217	!
2218	SUBROUTINE cpl2gath(champ_in, champ_out, klon, knon, iim, jjm, knindex)
2219
2220	! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
2221	! au coupleur.
2222	!
2223	!
2224	! input:
2225	! champ_in champ sur la grille gathere
2226	! knon nombre de points dans le domaine a traiter
2227	! knindex index des points de la surface a traiter
2228	! klon taille de la grille
2229	! iim,jjm dimension de la grille 2D
2230	!
2231	! output:
2232	! champ_out champ sur la grille 2D
2233	!
2234	! input
2235	integer :: klon, knon, iim, jjm
2236	real, dimension(iim,jjm+1) :: champ_in
2237	integer, dimension(klon) :: knindex
2238	! output
2239	real, dimension(klon) :: champ_out
2240	! local
2241	integer :: i, ig, j
2242	real, dimension(klon) :: tamp
2243	logical ,save :: check = .false.
2244
2245	ig = 1
2246	tamp(ig) = champ_in(1,1)
2247	do j = 2, jjm
2248	do i = 1, iim
2249	ig = ig + 1
2250	tamp(ig) = champ_in(i,j)
2251	enddo
2252	enddo
2253	ig = ig + 1
2254	tamp(ig) = champ_in(1,jjm+1)
2255
2256	do i = 1, knon
2257	ig = knindex(i)
2258	champ_out(i) = tamp(ig)
2259	enddo
2260
2261	END SUBROUTINE cpl2gath
2262	!
2263	!#########################################################################
2264	!
2265	SUBROUTINE albsno(klon, agesno,alb_neig_grid)
2266	IMPLICIT none
2267
2268	integer :: klon
2269	INTEGER, PARAMETER :: nvm = 8
2270	REAL, dimension(klon,nvm) :: veget
2271	REAL, DIMENSION(klon) :: alb_neig_grid, agesno
2272
2273	INTEGER :: i, nv
2274
2275	REAL, DIMENSION(nvm),SAVE :: init, decay
2276	REAL :: as
2277	DATA init /0.55, 0.14, 0.18, 0.29, 0.15, 0.15, 0.14, 0./
2278	DATA decay/0.30, 0.67, 0.63, 0.45, 0.40, 0.14, 0.06, 1./
2279
2280	veget = 0.
2281	veget(:,1) = 1. ! desert partout
2282	DO i = 1, klon
2283	alb_neig_grid(i) = 0.0
2284	ENDDO
2285	DO nv = 1, nvm
2286	DO i = 1, klon
2287	as = init(nv)+decay(nv)*EXP(-agesno(i)/5.)
2288	alb_neig_grid(i) = alb_neig_grid(i) + veget(i,nv)*as
2289	ENDDO
2290	ENDDO
2291
2292	END SUBROUTINE albsno
2293	!
2294	!#########################################################################
2295	!
2296
2297	SUBROUTINE fonte_neige( klon, knon, nisurf, dtime, &
2298	& tsurf, p1lay, cal, beta, coef1lay, ps, &
2299	& precip_rain, precip_snow, snow, qsol, &
2300	& radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
2301	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
2302	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
2303
2304	! Routine de traitement de la fonte de la neige dans le cas du traitement
2305	! de sol simplifié
2306	!
2307	! LF 03/2001
2308	! input:
2309	! knon nombre de points a traiter
2310	! nisurf surface a traiter
2311	! tsurf temperature de surface
2312	! p1lay pression 1er niveau (milieu de couche)
2313	! cal capacite calorifique du sol
2314	! beta evap reelle
2315	! coef1lay coefficient d'echange
2316	! ps pression au sol
2317	! precip_rain precipitations liquides
2318	! precip_snow precipitations solides
2319	! snow champs hauteur de neige
2320	! qsol humidite du sol
2321	! runoff runoff en cas de trop plein
2322	! petAcoef coeff. A de la resolution de la CL pour t
2323	! peqAcoef coeff. A de la resolution de la CL pour q
2324	! petBcoef coeff. B de la resolution de la CL pour t
2325	! peqBcoef coeff. B de la resolution de la CL pour q
2326	! radsol rayonnement net aus sol (LW + SW)
2327	! dif_grnd coeff. diffusion vers le sol profond
2328	!
2329	! output:
2330	! tsurf_new temperature au sol
2331	! fluxsens flux de chaleur sensible
2332	! fluxlat flux de chaleur latente
2333	! dflux_s derivee du flux de chaleur sensible / Ts
2334	! dflux_l derivee du flux de chaleur latente / Ts
2335	!
2336
2337	#include "YOETHF.inc"
2338	#include "FCTTRE.inc"
2339
2340	! Parametres d'entree
2341	integer, intent(IN) :: knon, nisurf, klon
2342	real , intent(IN) :: dtime
2343	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
2344	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
2345	real, dimension(klon), intent(IN) :: ps, q1lay
2346	real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
2347	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
2348	real, dimension(klon), intent(IN) :: radsol, dif_grnd
2349	real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
2350	real, dimension(klon), intent(INOUT) :: snow, qsol
2351
2352	! Parametres sorties
2353	real, dimension(klon), intent(INOUT):: tsurf_new, evap, fluxsens, fluxlat
2354	real, dimension(klon), intent(INOUT):: dflux_s, dflux_l
2355
2356	! Variables locales
2357	integer :: i
2358	real, dimension(klon) :: zx_mh, zx_nh, zx_oh
2359	real, dimension(klon) :: zx_mq, zx_nq, zx_oq
2360	real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
2361	real, dimension(klon) :: zx_sl, zx_k1
2362	real, dimension(klon) :: zx_q_0 , d_ts
2363	real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
2364	real :: bilan_f, fq_fonte
2365	REAL :: subli, fsno
2366	real, parameter :: t_grnd = 271.35, t_coup = 273.15
2367	!! PB temporaire en attendant mieux pour le modele de neige
2368	REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
2369	!
2370	logical, save :: check = .true.
2371	character (len = 20) :: modname = 'fonte_neige'
2372	logical, save :: neige_fond = .false.
2373	real, save :: max_eau_sol = 150.0
2374	character (len = 80) :: abort_message
2375	logical,save :: first = .true.,second=.false.
2376
2377	if (check) write(,)'Entree ', modname,' surface = ',nisurf
2378
2379	! Initialisations
2380	DO i = 1, knon
2381	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
2382	IF (thermcep) THEN
2383	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
2384	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
2385	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
2386	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
2387	zx_qs=MIN(0.5,zx_qs)
2388	zcor=1./(1.-retv*zx_qs)
2389	zx_qs=zx_qs*zcor
2390	zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
2391	& /RLVTT / zx_pkh(i)
2392	ELSE
2393	IF (tsurf(i).LT.t_coup) THEN
2394	zx_qs = qsats(tsurf(i)) / ps(i)
2395	zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
2396	& / zx_pkh(i)
2397	ELSE
2398	zx_qs = qsatl(tsurf(i)) / ps(i)
2399	zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
2400	& / zx_pkh(i)
2401	ENDIF
2402	ENDIF
2403	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
2404	zx_qsat(i) = zx_qs
2405	zx_coef(i) = coef1lay(i) &
2406	& * (1.0+SQRT(u1lay(i)2+v1lay(i)2)) &
2407	& * p1lay(i)/(RD*t1lay(i))
2408	ENDDO
2409
2410
2411	! === Calcul de la temperature de surface ===
2412	!
2413	! zx_sl = chaleur latente d'evaporation ou de sublimation
2414	!
2415	do i = 1, knon
2416	zx_sl(i) = RLVTT
2417	if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
2418	zx_k1(i) = zx_coef(i)
2419	enddo
2420
2421
2422	do i = 1, knon
2423	! Q
2424	zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
2425	zx_mq(i) = beta(i) * zx_k1(i) * &
2426	& (peqAcoef(i) - zx_qsat(i) &
2427	& + zx_dq_s_dt(i) * tsurf(i)) &
2428	& / zx_oq(i)
2429	zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
2430	& / zx_oq(i)
2431
2432	! H
2433	zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
2434	zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
2435	zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
2436	enddo
2437
2438	!
2439	! Y'a-t-il fonte de neige?
2440	!
2441	do i = 1, knon
2442	neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
2443	& .AND. tsurf_new(i) >= RTT)
2444	if (neige_fond) then
2445	tsurf_new(i) = RTT
2446	d_ts(i) = tsurf_new(i) - tsurf(i)
2447	! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
2448	! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
2449	!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
2450	!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
2451	evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
2452	fluxlat(i) = - evap(i) * zx_sl(i)
2453	fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
2454	! Derives des flux dF/dTs (W m-2 K-1):
2455	dflux_s(i) = zx_nh(i)
2456	dflux_l(i) = (zx_sl(i) * zx_nq(i))
2457	bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
2458	& dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
2459	& RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
2460	bilan_f = max(0., bilan_f)
2461	fq_fonte = bilan_f / zx_sl(i)
2462	snow(i) = max(0., snow(i) - fq_fonte * dtime)
2463	qsol(i) = qsol(i) + (fq_fonte * dtime)
2464	if (nisurf == is_ter) &
2465	& run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
2466	qsol(i) = min(qsol(i), max_eau_sol)
2467	endif
2468	enddo
2469
2470	END SUBROUTINE fonte_neige
2471	!
2472	!#########################################################################
2473	!
2474	END MODULE interface_surf

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

Download in other formats:

Original Format