Context Navigation

interface_surf.F90 @ 589

Last change on this file since 589 was 589, checked in by Laurent Fairhead, 19 years ago
Modifications pour le couplage carbone LOOP, PC LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 96.4 KB

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

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

Download in other formats:

Original Format