Context Navigation

interface_surf.F90 @ 757

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

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

Download in other formats:

Original Format