Context Navigation

interface_surf.F90 @ 1365

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

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

Download in other formats:

Original Format