Context Navigation

interface_surf.F90 @ 743

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

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

Download in other formats:

Original Format