Context Navigation

interface_surf.F90 @ 458

Last change on this file since 458 was 458, checked in by lmdzadmin, 21 years ago
Un peu de menage pour avoir moins de print dans les sorties LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 93.9 KB

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

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

Download in other formats:

Original Format