Context Navigation

interface_surf.F90 @ 457

Last change on this file since 457 was 457, checked in by lmdzadmin, 21 years ago
Correction bug tsurf_new=RTT - 1.8 apres fonte et passage routine soil Correction bug tsurf_new=RTT sic & lic IM/JLD/LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 92.5 KB

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

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

Download in other formats:

Original Format