Context Navigation

interface_surf.F90 @ 712

Last change on this file since 712 was 704, checked in by Laurent Fairhead, 19 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 122.8 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 n�essaire 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 adapt�!
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(knon), 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 parallel, only : pole_nord,pole_sud
1037	USE dimphy, klon_x=>klon
1038	IMPLICIT NONE
1039	! Cette routine sert d'interface entre le modele atmospherique et le
1040	! modele de sol continental. Appel a sechiba
1041	!
1042	! L. Fairhead 02/2000
1043	!
1044	! input:
1045	! itime numero du pas de temps
1046	! klon nombre total de points de grille
1047	! dtime pas de temps de la physique (en s)
1048	! nisurf index de la surface a traiter (1 = sol continental)
1049	! knon nombre de points de la surface a traiter
1050	! knindex index des points de la surface a traiter
1051	! rlon longitudes de la grille entiere
1052	! rlat latitudes de la grille entiere
1053	! pctsrf tableau des fractions de surface de chaque maille
1054	! debut logical: 1er appel a la physique (lire les restart)
1055	! lafin logical: dernier appel a la physique (ecrire les restart)
1056	! ok_veget logical: appel ou non au schema de surface continental
1057	! (si false calcul simplifie des fluxs sur les continents)
1058	! plev hauteur de la premiere couche (Pa)
1059	! u1_lay vitesse u 1ere couche
1060	! v1_lay vitesse v 1ere couche
1061	! temp_air temperature de l'air 1ere couche
1062	! spechum humidite specifique 1ere couche
1063	! epot_air temp pot de l'air
1064	! ccanopy concentration CO2 canopee
1065	! tq_cdrag cdrag
1066	! petAcoef coeff. A de la resolution de la CL pour t
1067	! peqAcoef coeff. A de la resolution de la CL pour q
1068	! petBcoef coeff. B de la resolution de la CL pour t
1069	! peqBcoef coeff. B de la resolution de la CL pour q
1070	! precip_rain precipitation liquide
1071	! precip_snow precipitation solide
1072	! lwdown flux IR descendant a la surface
1073	! swnet flux solaire net
1074	! swdown flux solaire entrant a la surface
1075	! tsurf temperature de surface
1076	! p1lay pression 1er niveau (milieu de couche)
1077	! ps pression au sol
1078	! radsol rayonnement net aus sol (LW + SW)
1079	!
1080	!
1081	! input/output
1082	! run_off ruissellement total
1083	!
1084	! output:
1085	! evap evaporation totale
1086	! fluxsens flux de chaleur sensible
1087	! fluxlat flux de chaleur latente
1088	! tsol_rad
1089	! tsurf_new temperature au sol
1090	! alb_new albedo
1091	! emis_new emissivite
1092	! z0_new surface roughness
1093	! qsurf air moisture at surface
1094
1095	! Parametres d'entree
1096	integer, intent(IN) :: itime
1097	integer, intent(IN) :: klon
1098	real, intent(IN) :: dtime
1099	real, intent(IN) :: date0
1100	integer, intent(IN) :: nisurf
1101	integer, intent(IN) :: knon
1102	integer, intent(IN) :: iim, jjm
1103	integer, dimension(klon), intent(IN) :: knindex
1104	logical, intent(IN) :: debut, lafin, ok_veget
1105	real, dimension(klon,nbsrf), intent(IN) :: pctsrf
1106	real, dimension(klon), intent(IN) :: rlon, rlat
1107	real, dimension(klon), intent(IN) :: cufi, cvfi
1108	real, dimension(klon), intent(IN) :: plev
1109	real, dimension(klon), intent(IN) :: u1_lay, v1_lay
1110	real, dimension(klon), intent(IN) :: temp_air, spechum
1111	real, dimension(klon), intent(IN) :: epot_air, ccanopy
1112	real, dimension(klon), intent(INOUT) :: tq_cdrag
1113	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
1114	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
1115	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
1116	real, dimension(klon), intent(IN) :: lwdown, swnet, swdown, ps
1117	!IM cf. JP +++
1118	real, dimension(klon) :: swdown_vrai
1119	!IM cf. JP ---
1120	real, dimension(klon), intent(IN) :: tsurf, p1lay
1121	real, dimension(klon), intent(IN) :: radsol
1122	! Parametres de sortie
1123	real, dimension(klon), intent(OUT):: evap, fluxsens, fluxlat, qsurf
1124	real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new, alblw
1125	real, dimension(klon), intent(OUT):: emis_new, z0_new
1126	real, dimension(klon), intent(OUT):: dflux_s, dflux_l
1127
1128	! Local
1129	!
1130	integer :: ii, ij, jj, igrid, ireal, i, index, iglob
1131	integer :: error
1132	character (len = 20) :: modname = 'interfsol'
1133	character (len = 80) :: abort_message
1134	logical,save :: check = .FALSE.
1135	real, dimension(klon) :: cal, beta, dif_grnd, capsol
1136	! type de couplage dans sechiba
1137	! character (len=10) :: coupling = 'implicit'
1138	! drapeaux controlant les appels dans SECHIBA
1139	! type(control_type), save :: control_in
1140	! Preserved albedo
1141	!IM cf. JP +++
1142	real, allocatable, dimension(:), save :: albedo_keep, zlev
1143	!IM cf. JP ---
1144	! coordonnees geographiques
1145	real, allocatable, dimension(:,:), save :: lalo
1146	! pts voisins
1147	integer,allocatable, dimension(:,:), save :: neighbours
1148	! fractions continents
1149	real,allocatable, dimension(:), save :: contfrac
1150	! resolution de la grille
1151	real, allocatable, dimension (:,:), save :: resolution
1152	! correspondance point n -> indices (i,j)
1153	integer, allocatable, dimension(:,:), save :: correspond
1154	! offset pour calculer les point voisins
1155	integer, dimension(8,3), save :: off_ini
1156	! Identifieurs des fichiers restart et histoire
1157	integer, save :: rest_id, hist_id
1158	integer, save :: rest_id_stom, hist_id_stom
1159	!
1160	real, allocatable, dimension (:,:), save :: lon_scat, lat_scat
1161
1162	logical, save :: lrestart_read = .true. , lrestart_write = .false.
1163
1164	real, dimension(klon):: snow
1165	real, dimension(knon,2) :: albedo_out
1166	! Pb de nomenclature
1167	real, dimension(klon) :: petA_orc, peqA_orc
1168	real, dimension(klon) :: petB_orc, peqB_orc
1169	! Pb de correspondances de grilles
1170	integer, dimension(:), save, allocatable :: ig, jg
1171	integer :: indi, indj
1172	integer, save, allocatable,dimension(:) :: ktindex
1173	REAL, dimension(klon) :: bidule
1174	! Essai cdrag
1175	real, dimension(klon) :: cdrag
1176	integer :: jjb,jje,ijb,ije
1177	INTEGER,SAVE :: offset
1178	REAL, dimension(klon2) :: rlon_g,rlat_g
1179	INTEGER, SAVE :: orch_comm
1180	#include "temps.inc"
1181	#include "YOMCST.inc"
1182	#include "iniprint.h"
1183
1184	if (check) write(lunout,*)'Entree ', modname
1185	if (check) write(lunout,*)'ok_veget = ',ok_veget
1186
1187
1188
1189	! initialisation
1190
1191	if (debut) then
1192	ALLOCATE(ktindex(klon))
1193	IF ( .NOT. allocated(albedo_keep)) THEN
1194	ALLOCATE(albedo_keep(klon))
1195	ALLOCATE(zlev(klon))
1196	ENDIF
1197	! Pb de correspondances de grilles
1198	allocate(ig(klon))
1199	allocate(jg(klon))
1200	ig(1) = 1
1201	jg(1) = 1
1202	indi = 0
1203	indj = 2
1204	do igrid = 2, klon - 1
1205	indi = indi + 1
1206	if ( indi > iim) then
1207	indi = 1
1208	indj = indj + 1
1209	endif
1210	ig(igrid) = indi
1211	jg(igrid) = indj
1212	enddo
1213	ig(klon) = 1
1214	jg(klon) = jjm + 1
1215
1216	if ((.not. allocated(lalo))) then
1217	allocate(lalo(knon,2), stat = error)
1218	if (error /= 0) then
1219	abort_message='Pb allocation lalo'
1220	call abort_gcm(modname,abort_message,1)
1221	endif
1222	endif
1223	if ((.not. allocated(lon_scat))) then
1224	allocate(lon_scat(iim,jjm+1), stat = error)
1225	if (error /= 0) then
1226	abort_message='Pb allocation lon_scat'
1227	call abort_gcm(modname,abort_message,1)
1228	endif
1229	endif
1230	if ((.not. allocated(lat_scat))) then
1231	allocate(lat_scat(iim,jjm+1), stat = error)
1232	if (error /= 0) then
1233	abort_message='Pb allocation lat_scat'
1234	call abort_gcm(modname,abort_message,1)
1235	endif
1236	endif
1237	lon_scat = 0.
1238	lat_scat = 0.
1239	do igrid = 1, knon
1240	index = knindex(igrid)
1241	lalo(igrid,2) = rlon(index)
1242	lalo(igrid,1) = rlat(index)
1243
1244	enddo
1245
1246
1247
1248	Call GatherField(rlon,rlon_g,1)
1249	Call GatherField(rlat,rlat_g,1)
1250
1251	IF (phy_rank==0) THEN
1252	index = 1
1253	do jj = 2, jjm
1254	do ij = 1, iim
1255	index = index + 1
1256	lon_scat(ij,jj) = rlon_g(index)
1257	lat_scat(ij,jj) = rlat_g(index)
1258	enddo
1259	enddo
1260	lon_scat(:,1) = lon_scat(:,2)
1261	lat_scat(:,1) = rlat_g(1)
1262	lon_scat(:,jjm+1) = lon_scat(:,2)
1263	lat_scat(:,jjm+1) = rlat_g(klon2)
1264	ENDIF
1265
1266
1267	!
1268	! Allouer et initialiser le tableau des voisins et des fraction de continents
1269	!
1270	if ( (.not.allocated(neighbours))) THEN
1271	allocate(neighbours(knon,8), stat = error)
1272	if (error /= 0) then
1273	abort_message='Pb allocation neighbours'
1274	call abort_gcm(modname,abort_message,1)
1275	endif
1276	endif
1277	neighbours = -1.
1278	if (( .not. allocated(contfrac))) then
1279	allocate(contfrac(knon), stat = error)
1280	if (error /= 0) then
1281	abort_message='Pb allocation contfrac'
1282	call abort_gcm(modname,abort_message,1)
1283	endif
1284	endif
1285
1286	do igrid = 1, knon
1287	ireal = knindex(igrid)
1288	contfrac(igrid) = pctsrf(ireal,is_ter)
1289	enddo
1290
1291
1292	CALL Init_neighbours(iim,jjm,knon,neighbours,knindex,pctsrf(:,is_ter))
1293
1294	!
1295	! Allocation et calcul resolutions
1296	IF ( (.NOT.ALLOCATED(resolution))) THEN
1297	ALLOCATE(resolution(knon,2), stat = error)
1298	if (error /= 0) then
1299	abort_message='Pb allocation resolution'
1300	call abort_gcm(modname,abort_message,1)
1301	endif
1302	ENDIF
1303	do igrid = 1, knon
1304	ij = knindex(igrid)
1305	resolution(igrid,1) = cufi(ij)
1306	resolution(igrid,2) = cvfi(ij)
1307	enddo
1308
1309	endif ! (fin debut)
1310
1311	!
1312	! Appel a la routine sols continentaux
1313	!
1314	if (lafin) lrestart_write = .true.
1315	if (check) write(lunout,*)'lafin ',lafin,lrestart_write
1316
1317	petA_orc = petBcoef * dtime
1318	petB_orc = petAcoef
1319	peqA_orc = peqBcoef * dtime
1320	peqB_orc = peqAcoef
1321
1322	cdrag = 0.
1323	cdrag(1:knon) = tq_cdrag(1:knon)
1324
1325	!IM cf. JP +++
1326	! zlev(1:knon) = (100.plev(1:knon))/((ps(1:knon)/287.05temp_air(1:knon))*9.80665)
1327	zlev(1:knon) = (100.plev(1:knon))/((ps(1:knon)/RDtemp_air(1:knon))*RG)
1328	!IM cf. JP ---
1329
1330
1331	! PF et PASB
1332	! where(cdrag > 0.01)
1333	! cdrag = 0.01
1334	! endwhere
1335	! write(,)'Cdrag = ',minval(cdrag),maxval(cdrag)
1336
1337	!
1338	! Init Orchidee
1339	!
1340	! if (pole_nord) then
1341	! offset=0
1342	! ktindex(:)=ktindex(:)+iim-1
1343	! else
1344	! offset = klon_begin-1+iim-1
1345	! ktindex(:)=ktindex(:)+MOD(offset,iim)
1346	! offset=offset-MOD(offset,iim)
1347	! endif
1348
1349	PRINT *,'ORCHIDEE ------> KNON : ',knon
1350
1351
1352	if (debut) then
1353	CALL Get_orchidee_communicator(knon,orch_comm)
1354	IF (knon /=0) THEN
1355	CALL Init_orchidee_index(iim,knon,orch_comm,knindex,offset,ktindex)
1356
1357	call intersurf_main (itime+itau_phy-1, iim, jjm+1,offset, knon, ktindex, &
1358	& orch_comm,dtime, lrestart_read, lrestart_write, lalo, &
1359	& contfrac, neighbours, resolution, date0, &
1360	& zlev, u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
1361	& cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
1362	& precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
1363	& evap, fluxsens, fluxlat, coastalflow, riverflow, &
1364	& tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
1365	& lon_scat, lat_scat)
1366
1367	ENDIF
1368	!IM cf. JP +++
1369	albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
1370	!IM cf. JP ---
1371
1372	endif
1373
1374	!IM cf. JP +++
1375	! swdown_vrai(1:knon) = swnet(1:knon)/(1. - albedo_keep(1:knon))
1376	swdown_vrai(1:knon) = swdown(1:knon)
1377
1378	!IM cf. JP ---
1379	IF (knon /=0) THEN
1380
1381	call intersurf_main (itime+itau_phy, iim, jjm+1,offset, knon, ktindex, &
1382	& orch_comm,dtime, lrestart_read, lrestart_write, lalo, &
1383	& contfrac, neighbours, resolution, date0, &
1384	& zlev, u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
1385	& cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
1386	& precip_rain, precip_snow, lwdown, swnet, swdown_vrai, ps, &
1387	& evap, fluxsens, fluxlat, coastalflow, riverflow, &
1388	& tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
1389	& lon_scat, lat_scat)
1390
1391	ENDIF
1392	!IM cf. JP +++
1393	albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
1394	!IM cf. JP ---
1395
1396	bidule=0.
1397	bidule(1:knon)=riverflow(1:knon)
1398	call gath2cpl(bidule, tmp_rriv, klon, knon,iim,jjphy_nb,knindex)
1399	bidule=0.
1400	bidule(1:knon)=coastalflow(1:knon)
1401	call gath2cpl(bidule, tmp_rcoa, klon, knon,iim,jjphy_nb,knindex)
1402	alb_new(1:knon) = albedo_out(1:knon,1)
1403	alblw(1:knon) = albedo_out(1:knon,2)
1404
1405
1406	! Convention orchidee: positif vers le haut
1407	fluxsens(1:knon) = -1. * fluxsens(1:knon)
1408	fluxlat(1:knon) = -1. * fluxlat(1:knon)
1409
1410	! evap = -1. * evap
1411
1412	if (debut) lrestart_read = .false.
1413
1414	END SUBROUTINE interfsol
1415
1416	SUBROUTINE Init_orchidee_index(iim,knon,orch_comm,knindex,offset,ktindex)
1417	USE dimphy
1418	IMPLICIT NONE
1419	INTEGER,INTENT(IN) :: iim
1420	INTEGER,INTENT(IN) :: knon
1421	INTEGER,INTENT(IN) :: orch_comm
1422	INTEGER,INTENT(IN) :: knindex(knon)
1423	INTEGER,INTENT(OUT) :: offset
1424	INTEGER,INTENT(OUT) :: ktindex(knon)
1425
1426	#ifdef CPP_PARA
1427	INCLUDE 'mpif.h'
1428	INTEGER :: status(MPI_STATUS_SIZE)
1429	#endif
1430	INTEGER :: MyLastPoint
1431	INTEGER :: LastPoint
1432	INTEGER :: mpi_rank
1433	INTEGER :: mpi_size
1434	INTEGER :: ierr
1435
1436	MyLastPoint=klon_begin-1+knindex(knon)+iim-1
1437
1438	IF (.NOT. monocpu) THEN
1439	#ifdef CPP_PARA
1440	call MPI_COMM_SIZE(orch_comm,mpi_size,ierr)
1441	call MPI_COMM_RANK(orch_comm,mpi_rank,ierr)
1442	#endif
1443	ELSE
1444	mpi_rank=0
1445	mpi_size=1
1446	ENDIF
1447
1448	IF (.NOT. monocpu) THEN
1449	IF (mpi_rank /= 0) then
1450	#ifdef CPP_PARA
1451	CALL MPI_RECV(LastPoint,1,MPI_INTEGER,mpi_rank-1,1234,orch_comm,status,ierr)
1452	#endif
1453	ENDIF
1454
1455	IF (mpi_rank /= mpi_size-1) THEN
1456	#ifdef CPP_PARA
1457	CALL MPI_SEND(MyLastPoint,1,MPI_INTEGER,mpi_rank+1,1234,orch_comm,ierr)
1458	#endif
1459	ENDIF
1460	ENDIF
1461
1462	IF (mpi_rank==0) THEN
1463	offset=0
1464	ELSE
1465	offset=LastPoint-MOD(LastPoint,iim)
1466	ENDIF
1467
1468	ktindex(:)=knindex(:)+(klon_begin+iim-1)-offset-1
1469
1470
1471	END SUBROUTINE Init_orchidee_index
1472
1473
1474	SUBROUTINE Get_orchidee_communicator(knon,orch_comm)
1475	USE dimphy, only : phy_rank
1476	USE parallel, only : COMM_LMDZ
1477	IMPLICIT NONE
1478	#ifdef CPP_PARA
1479	include 'mpif.h'
1480	#endif
1481	INTEGER,INTENT(IN) :: knon
1482	INTEGER,INTENT(OUT) :: orch_comm
1483
1484	INTEGER :: color
1485	INTEGER :: ierr
1486
1487	IF (knon==0) THEN
1488	color = 0
1489	ELSE
1490	color = 1
1491	ENDIF
1492
1493	#ifdef CPP_PARA
1494	CALL MPI_COMM_SPLIT(COMM_LMDZ,color,phy_rank,orch_comm,ierr)
1495	#endif
1496
1497	END SUBROUTINE Get_orchidee_communicator
1498
1499
1500	SUBROUTINE Init_neighbours(iim,jjm,knon,neighbours,ktindex,pctsrf)
1501	USE parallel,only : COMM_LMDZ
1502	USE dimphy
1503	IMPLICIT NONE
1504	#ifdef CPP_PARA
1505	include 'mpif.h'
1506	#endif
1507	INTEGER :: iim,jjm
1508	INTEGER :: knon
1509	INTEGER :: neighbours(knon,8)
1510	INTEGER :: ktindex(knon)
1511	REAL :: pctsrf(klon)
1512
1513	INTEGER :: knon_nb(0:phy_size-1)
1514	INTEGER,DIMENSION(0:phy_size-1) :: displs,sendcount
1515	INTEGER,ALLOCATABLE :: ktindex_g(:)
1516	REAL*8 :: pctsrf_g(klon2)
1517	INTEGER,ALLOCATABLE ::neighbours_g(:,:)
1518	INTEGER :: knon_g
1519	REAL*8 :: correspond(iim,jjm+1)
1520	INTEGER :: i,igrid,jj,ij,iglob,ierr,ireal,index
1521	integer, dimension(8,3) :: off_ini
1522	integer, dimension(8) :: offset
1523	INTEGER :: ktindex_p(knon)
1524
1525	IF (monocpu) THEN
1526	knon_nb(:)=knon
1527	ELSE
1528
1529	#ifdef CPP_PARA
1530	CALL MPI_GATHER(knon,1,MPI_INTEGER,knon_nb,1,MPI_INTEGER,0,COMM_LMDZ,ierr)
1531	#endif
1532
1533	ENDIF
1534
1535	IF (phy_rank==0) THEN
1536	knon_g=sum(knon_nb(:))
1537	ALLOCATE(ktindex_g(knon_g))
1538	ALLOCATE(neighbours_g(knon_g,8))
1539	neighbours_g(:,:)=-1
1540	displs(0)=0
1541	DO i=1,phy_size-1
1542	displs(i)=displs(i-1)+knon_nb(i-1)
1543	ENDDO
1544	ENDIF
1545
1546	ktindex_p(:)=ktindex(:)+klon_begin-1+iim-1
1547
1548	IF (monocpu) THEN
1549	ktindex_g(:)=ktindex_p(:)
1550	ELSE
1551
1552	#ifdef CPP_PARA
1553	CALL MPI_GATHERV(ktindex_p,knon,MPI_INTEGER,ktindex_g,knon_nb,displs,MPI_INTEGER,0,COMM_LMDZ,ierr)
1554	#endif
1555
1556	ENDIF
1557
1558	CALL GatherField(pctsrf,pctsrf_g,1)
1559
1560	IF (phy_rank==0) THEN
1561	! Initialisation des offset
1562	!
1563	! offset bord ouest
1564	off_ini(1,1) = - iim ; off_ini(2,1) = - iim + 1; off_ini(3,1) = 1
1565	off_ini(4,1) = iim + 1; off_ini(5,1) = iim ; off_ini(6,1) = 2 * iim - 1
1566	off_ini(7,1) = iim -1 ; off_ini(8,1) = - 1
1567	! offset point normal
1568	off_ini(1,2) = - iim ; off_ini(2,2) = - iim + 1; off_ini(3,2) = 1
1569	off_ini(4,2) = iim + 1; off_ini(5,2) = iim ; off_ini(6,2) = iim - 1
1570	off_ini(7,2) = -1 ; off_ini(8,2) = - iim - 1
1571	! offset bord est
1572	off_ini(1,3) = - iim; off_ini(2,3) = - 2 * iim + 1; off_ini(3,3) = - iim + 1
1573	off_ini(4,3) = 1 ; off_ini(5,3) = iim ; off_ini(6,3) = iim - 1
1574	off_ini(7,3) = -1 ; off_ini(8,3) = - iim - 1
1575	!
1576	!
1577	! Attention aux poles
1578	!
1579	do igrid = 1, knon_g
1580	index = ktindex_g(igrid)
1581	jj = int((index - 1)/iim) + 1
1582	ij = index - (jj - 1) * iim
1583	correspond(ij,jj) = igrid
1584	enddo
1585
1586	do igrid = 1, knon_g
1587	iglob = ktindex_g(igrid)
1588	if (mod(iglob, iim) == 1) then
1589	offset = off_ini(:,1)
1590	else if(mod(iglob, iim) == 0) then
1591	offset = off_ini(:,3)
1592	else
1593	offset = off_ini(:,2)
1594	endif
1595	do i = 1, 8
1596	index = iglob + offset(i)
1597	ireal = (min(max(1, index - iim + 1), klon2))
1598	if (pctsrf_g(ireal) > EPSFRA) then
1599	jj = int((index - 1)/iim) + 1
1600	ij = index - (jj - 1) * iim
1601	neighbours_g(igrid, i) = correspond(ij, jj)
1602	endif
1603	enddo
1604	enddo
1605
1606	! DO i=0,phy_size-1
1607	! displs(i)=displs(i)*8
1608	! sendcount(i)=knon_nb(i)*8
1609	! ENDDO
1610
1611	ENDIF
1612
1613	DO i=1,8
1614	IF (monocpu) THEN
1615	neighbours(:,i)=neighbours_g(:,i)
1616	ELSE
1617	#ifdef CPP_PARA
1618	CALL MPI_SCATTERV(neighbours_g(:,i),knon_nb,displs,MPI_INTEGER,neighbours(:,i),knon,MPI_INTEGER,0,COMM_LMDZ,ierr)
1619	#endif
1620	ENDIF
1621	ENDDO
1622
1623	END SUBROUTINE Init_neighbours
1624	#endif
1625	!
1626	!#########################################################################
1627	!
1628	SUBROUTINE interfoce_cpl(itime, dtime, cumul, &
1629	& klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
1630	& ocean, npas, nexca, debut, lafin, &
1631	& swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
1632	& fluxlat, fluxsens, fder, albsol, taux, tauy, &
1633	! -- LOOP
1634	& windsp, &
1635	! -- LOOP
1636	& zmasq, &
1637	& tsurf_new, alb_new, &
1638	& pctsrf_new)
1639
1640	USE ioipsl
1641	USE dimphy, only : jjphy_nb, iiphy_begin,iiphy_end,phy_rank,phy_size, monocpu
1642	USE iophy
1643	#ifdef CPP_PARA
1644	USE parallel, only: pole_nord,pole_sud,COMM_LMDZ
1645	#endif
1646	#ifdef CPP_PSMILE
1647	USE oasis
1648	#endif
1649	USE write_field_phy
1650	implicit none
1651	#include "indicesol.inc"
1652	#include "YOMCST.inc"
1653	! Cette routine sert d'interface entre le modele atmospherique et un
1654	! coupleur avec un modele d'ocean 'complet' derriere
1655	!
1656	! Le modele de glace qu'il est prevu d'utiliser etant couple directement a
1657	! l'ocean presentement, on va passer deux fois dans cette routine par pas de
1658	! temps physique, une fois avec les points oceans et l'autre avec les points
1659	! glace. A chaque pas de temps de couplage, la lecture des champs provenant
1660	! du coupleur se fera "dans" l'ocean et l'ecriture des champs a envoyer
1661	! au coupleur "dans" la glace. Il faut donc des tableaux de travail "tampons"
1662	! dimensionnes sur toute la grille qui remplissent les champs sur les
1663	! domaines ocean/glace quand il le faut. Il est aussi necessaire que l'index
1664	! ocean soit traiter avant l'index glace (sinon tout intervertir)
1665	!
1666	!
1667	! L. Fairhead 02/2000
1668	!
1669	! input:
1670	! itime numero du pas de temps
1671	! iim, jjm nbres de pts de grille
1672	! dtime pas de temps de la physique
1673	! klon nombre total de points de grille
1674	! nisurf index de la surface a traiter (1 = sol continental)
1675	! pctsrf tableau des fractions de surface de chaque maille
1676	! knon nombre de points de la surface a traiter
1677	! knindex index des points de la surface a traiter
1678	! rlon longitudes
1679	! rlat latitudes
1680	! debut logical: 1er appel a la physique
1681	! lafin logical: dernier appel a la physique
1682	! ocean type d'ocean
1683	! nexca frequence de couplage
1684	! swdown flux solaire entrant a la surface
1685	! lwdown flux IR net a la surface
1686	! precip_rain precipitation liquide
1687	! precip_snow precipitation solide
1688	! evap evaporation
1689	! tsurf temperature de surface
1690	! fder derivee dF/dT
1691	! albsol albedo du sol (coherent avec swdown)
1692	! taux tension de vent en x
1693	! tauy tension de vent en y
1694	! -- LOOP
1695	! windsp module du vent a 10m
1696	! -- LOOP
1697	! nexca frequence de couplage
1698	! zmasq masque terre/ocean
1699	!
1700	!
1701	! output:
1702	! tsurf_new temperature au sol
1703	! alb_new albedo
1704	! pctsrf_new nouvelle repartition des surfaces
1705	! alb_ice albedo de la glace
1706	!
1707
1708
1709	! Parametres d'entree
1710	integer, intent(IN) :: itime
1711	integer, intent(IN) :: iim, jjm
1712	real, intent(IN) :: dtime
1713	integer, intent(IN) :: klon
1714	integer, intent(IN) :: nisurf
1715	integer, intent(IN) :: knon
1716	real, dimension(klon,nbsrf), intent(IN) :: pctsrf
1717	integer, dimension(klon), intent(in) :: knindex
1718	logical, intent(IN) :: debut, lafin
1719	real, dimension(klon), intent(IN) :: rlon, rlat
1720	character (len = 6) :: ocean
1721	real, dimension(klon), intent(IN) :: lwdown, swdown
1722	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
1723	real, dimension(klon), intent(IN) :: tsurf, fder, albsol, taux, tauy
1724	! -- LOOP
1725	real, dimension(klon), intent(IN) :: windsp
1726	! -- LOOP
1727	INTEGER :: nexca, npas, kstep
1728	real, dimension(klon), intent(IN) :: zmasq
1729	real, dimension(klon), intent(IN) :: fluxlat, fluxsens
1730	logical, intent(IN) :: cumul
1731	real, dimension(klon), intent(INOUT) :: evap
1732
1733	! Parametres de sortie
1734	real, dimension(klon), intent(OUT):: tsurf_new, alb_new
1735	real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
1736
1737	! Variables locales
1738	integer :: j, error, sum_error, ig, cpl_index,i
1739	! -- LOOP
1740	INTEGER :: nsrf
1741	! -- LOOP
1742	character (len = 20) :: modname = 'interfoce_cpl'
1743	character (len = 80) :: abort_message
1744	logical,save :: check = .FALSE.
1745	! variables pour moyenner les variables de couplage
1746	real, allocatable, dimension(:,:),save :: cpl_sols, cpl_nsol, cpl_rain
1747	real, allocatable, dimension(:,:),save :: cpl_snow, cpl_evap, cpl_tsol
1748	real, allocatable, dimension(:,:),save :: cpl_fder, cpl_albe, cpl_taux
1749	! -- LOOP
1750	real, allocatable, dimension(:,:),save :: cpl_windsp
1751	! -- LOOP
1752	real, allocatable, dimension(:,:),save :: cpl_tauy
1753	REAL, ALLOCATABLE, DIMENSION(:,:),SAVE :: cpl_rriv, cpl_rcoa, cpl_rlic
1754	!!$
1755	! variables tampons avant le passage au coupleur
1756	real, allocatable, dimension(:,:,:),save :: tmp_sols, tmp_nsol, tmp_rain
1757	real, allocatable, dimension(:,:,:),save :: tmp_snow, tmp_evap, tmp_tsol
1758	real, allocatable, dimension(:,:,:),save :: tmp_fder, tmp_albe, tmp_taux
1759	! -- LOOP
1760	real, allocatable, dimension(:,:,:),save :: tmp_windsp
1761	! -- LOOP
1762	!!$ real, allocatable, dimension(:,:,:),save :: tmp_tauy, tmp_rriv, tmp_rcoa
1763	REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: tmp_tauy
1764	! variables a passer au coupleur
1765	!ym real, dimension(iim, jjm+1) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice
1766	!ym real, dimension(iim, jjm+1) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
1767	!ym REAL, DIMENSION(iim, jjm+1) :: wri_evap_sea, wri_rcoa, wri_rriv
1768	!ym REAL, DIMENSION(iim, jjm+1) :: wri_rain, wri_snow, wri_taux, wri_tauy
1769	!ym REAL, DIMENSION(iim, jjm+1) :: wri_calv
1770	!ym REAL, DIMENSION(iim, jjm+1) :: wri_tauxx, wri_tauyy, wri_tauzz
1771	!ym REAL, DIMENSION(iim, jjm+1) :: tmp_lon, tmp_lat
1772
1773	real, dimension(iim, jjphy_nb) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice
1774	real, dimension(iim, jjphy_nb) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
1775	REAL, DIMENSION(iim, jjphy_nb) :: wri_evap_sea, wri_rcoa, wri_rriv
1776	REAL, DIMENSION(iim, jjphy_nb) :: wri_rain, wri_snow, wri_taux, wri_tauy
1777	REAL, DIMENSION(iim, jjphy_nb) :: wri_calv
1778	REAL, DIMENSION(iim, jjphy_nb) :: wri_tauxx, wri_tauyy, wri_tauzz
1779	REAL, DIMENSION(iim, jjphy_nb) :: tmp_lon, tmp_lat
1780	REAL, DIMENSION(iim, jjphy_nb) :: wri_windsp
1781
1782	! variables relues par le coupleur
1783	! read_sic = fraction de glace
1784	! read_sit = temperature de glace
1785	real, allocatable, dimension(:,:),save :: read_sst, read_sic, read_sit
1786	real, allocatable, dimension(:,:),save :: read_alb_sic
1787	! variable tampon
1788	real, dimension(klon) :: tamp_sic
1789	! sauvegarde des fractions de surface d'un pas de temps a l'autre apres
1790	! l'avoir lu
1791	real, allocatable,dimension(:,:),save :: pctsrf_sav
1792	real, dimension(iim, jjphy_nb, 3) :: tamp_srf
1793	integer, allocatable, dimension(:), save :: tamp_ind
1794	real, allocatable, dimension(:,:),save :: tamp_zmasq
1795	real, dimension(iim, jjphy_nb) :: deno
1796	integer :: idtime
1797	integer, allocatable,dimension(:),save :: unity
1798	!
1799	logical, save :: first_appel = .true.
1800	logical,save :: print
1801	!maf
1802	! variables pour avoir une sortie IOIPSL des champs echanges
1803	CHARACTER*80,SAVE :: clintocplnam, clfromcplnam
1804	INTEGER, SAVE :: jf,nhoridct,nidct
1805	INTEGER, SAVE :: nhoridcs,nidcs
1806	INTEGER :: ndexct(iim(jjm+1)),ndexcs(iim(jjm+1))
1807	REAL :: zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
1808	INTEGER,save :: idayref
1809	!med integer :: itau_w
1810	integer,save :: itau_w
1811	! -- LOOP
1812	integer :: nb_interf_cpl
1813	! -- LOOP
1814
1815	real :: Up,Down
1816	integer :: ierr
1817	integer :: il_time_secs
1818	real :: tmp_field(klon)
1819
1820	#include "param_cou.h"
1821	#include "inc_cpl.h"
1822	#include "temps.inc"
1823	#include "iniprint.h"
1824
1825	#ifdef CPP_PARA
1826	include 'mpif.h'
1827	integer :: status(MPI_STATUS_SIZE)
1828	#endif
1829
1830	!
1831	! Initialisation
1832	!
1833	if (check) write(,)'Entree ',modname,'nisurf = ',nisurf
1834
1835	if (first_appel) then
1836	error = 0
1837	allocate(unity(klon), stat = error)
1838	if ( error /=0) then
1839	abort_message='Pb allocation variable unity'
1840	call abort_gcm(modname,abort_message,1)
1841	endif
1842	allocate(pctsrf_sav(klon,nbsrf), stat = error)
1843	if ( error /=0) then
1844	abort_message='Pb allocation variable pctsrf_sav'
1845	call abort_gcm(modname,abort_message,1)
1846	endif
1847	pctsrf_sav = 0.
1848
1849	do ig = 1, klon
1850	unity(ig) = ig
1851	enddo
1852	sum_error = 0
1853	allocate(cpl_sols(klon,2), stat = error); sum_error = sum_error + error
1854	allocate(cpl_nsol(klon,2), stat = error); sum_error = sum_error + error
1855	allocate(cpl_rain(klon,2), stat = error); sum_error = sum_error + error
1856	allocate(cpl_snow(klon,2), stat = error); sum_error = sum_error + error
1857	allocate(cpl_evap(klon,2), stat = error); sum_error = sum_error + error
1858	allocate(cpl_tsol(klon,2), stat = error); sum_error = sum_error + error
1859	allocate(cpl_fder(klon,2), stat = error); sum_error = sum_error + error
1860	allocate(cpl_albe(klon,2), stat = error); sum_error = sum_error + error
1861	allocate(cpl_taux(klon,2), stat = error); sum_error = sum_error + error
1862	! -- LOOP
1863	allocate(cpl_windsp(klon,2), stat = error); sum_error = sum_error + error
1864	! -- LOOP
1865	allocate(cpl_tauy(klon,2), stat = error); sum_error = sum_error + error
1866	!ym ALLOCATE(cpl_rriv(iim,jjm+1), stat=error); sum_error = sum_error + error
1867	!ym ALLOCATE(cpl_rcoa(iim,jjm+1), stat=error); sum_error = sum_error + error
1868	!ym ALLOCATE(cpl_rlic(iim,jjm+1), stat=error); sum_error = sum_error + error
1869	ALLOCATE(cpl_rriv(iim,jjphy_nb), stat=error); sum_error = sum_error + error
1870	ALLOCATE(cpl_rcoa(iim,jjphy_nb), stat=error); sum_error = sum_error + error
1871	ALLOCATE(cpl_rlic(iim,jjphy_nb), stat=error); sum_error = sum_error + error
1872
1873
1874	!!
1875	!ym allocate(read_sst(iim, jjm+1), stat = error); sum_error = sum_error + error
1876	!ym allocate(read_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
1877	!ym allocate(read_sit(iim, jjm+1), stat = error); sum_error = sum_error + error
1878	!ym allocate(read_alb_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
1879	allocate(read_sst(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1880	allocate(read_sic(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1881	allocate(read_sit(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1882	allocate(read_alb_sic(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1883	read_sst=0.
1884	read_sic=0.
1885	read_sit=0.
1886	read_alb_sic=0.
1887	if (sum_error /= 0) then
1888	abort_message='Pb allocation variables couplees'
1889	call abort_gcm(modname,abort_message,1)
1890	endif
1891	cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
1892	cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
1893	cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.; cpl_rlic = 0.
1894	! -- LOOP
1895	cpl_windsp = 0.
1896	! -- LOOP
1897
1898	sum_error = 0
1899	allocate(tamp_ind(klon), stat = error); sum_error = sum_error + error
1900	!ym allocate(tamp_zmasq(iim, jjm+1), stat = error); sum_error = sum_error + error
1901	allocate(tamp_zmasq(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1902	tamp_zmasq=1.
1903
1904	do ig = 1, klon
1905	tamp_ind(ig) = ig
1906	enddo
1907	call gath2cpl(zmasq, tamp_zmasq, klon, klon, iim, jjphy_nb, tamp_ind)
1908	!
1909	! initialisation couplage
1910	!
1911	idtime = int(dtime)
1912	#ifdef CPP_COUPLE
1913	#ifdef CPP_PSMILE
1914	CALL inicma(iim, (jjm+1))
1915	#else
1916	if (.not. monocpu) then
1917	abort_message='coupleur parallele uniquement avec PSMILE'
1918	call abort_gcm(modname,abort_message,1)
1919	endif
1920	call inicma(npas , nexca, idtime,(jjm+1)*iim)
1921	#endif
1922	#endif
1923	!
1924	! initialisation sorties netcdf
1925	!
1926	!ym IO de check deconnect�pour le moment en //
1927	IF (monocpu) THEN
1928	idayref = day_ini
1929	CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
1930	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
1931	DO i = 1, iim
1932	zx_lon(i,1) = rlon(i+1)
1933	zx_lon(i,jjm+1) = rlon(i+1)
1934	ENDDO
1935	CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
1936	clintocplnam="cpl_atm_tauflx"
1937	CALL histbeg(clintocplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
1938	& itau_phy,zjulian,dtime,nhoridct,nidct)
1939	! no vertical axis
1940	CALL histdef(nidct, 'tauxe','tauxe', &
1941	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1942	CALL histdef(nidct, 'tauyn','tauyn', &
1943	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1944	CALL histdef(nidct, 'tmp_lon','tmp_lon', &
1945	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1946	CALL histdef(nidct, 'tmp_lat','tmp_lat', &
1947	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1948	DO jf=1,jpflda2o1 + jpflda2o2
1949	CALL histdef(nidct, cl_writ(jf),cl_writ(jf), &
1950	& "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1951	END DO
1952	CALL histend(nidct)
1953	CALL histsync(nidct)
1954
1955	clfromcplnam="cpl_atm_sst"
1956	CALL histbeg(clfromcplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
1957	& 0,zjulian,dtime,nhoridcs,nidcs)
1958	! no vertical axis
1959	DO jf=1,jpfldo2a
1960	CALL histdef(nidcs, cl_read(jf),cl_read(jf), &
1961	& "-",iim, jjm+1, nhoridcs, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1962	END DO
1963	CALL histend(nidcs)
1964	CALL histsync(nidcs)
1965
1966	ENDIF ! monocpu
1967
1968	! pour simuler la fonte des glaciers antarctiques
1969	!
1970	!ym => pour le moment, c'est en commentaire, donc je squizze
1971
1972	!ym surf_maille = (4. * rpi * ra*2) / (iim (jjm +1))
1973	!ym ALLOCATE(coeff_iceberg(iim,jjm+1), stat=error)
1974	!ym if (error /= 0) then
1975	!ym abort_message='Pb allocation variable coeff_iceberg'
1976	!ym call abort_gcm(modname,abort_message,1)
1977	!ym endif
1978	!ym open (12,file='flux_iceberg',form='formatted',status='old')
1979	!ym read (12,*) coeff_iceberg
1980	!ym close (12)
1981	!ym num_antarctic = max(1, count(coeff_iceberg > 0))
1982
1983	first_appel = .false.
1984	endif ! fin if (first_appel)
1985
1986	! Initialisations
1987
1988	! calcul des fluxs a passer
1989	! -- LOOP
1990	nb_interf_cpl = nb_interf_cpl + 1
1991	if (check) write(lunout,*)'passage dans interface_surf.F90 : ',nb_interf_cpl
1992	! -- LOOP
1993	cpl_index = 1
1994	if (nisurf == is_sic) cpl_index = 2
1995	if (cumul) then
1996	! -- LOOP
1997	if (check) write(lunout,*)'passage dans cumul '
1998	if (check) write(lunout,*)'valeur de cpl_index ', cpl_index
1999	! -- LOOP
2000	if (check) write(,) modname, 'cumul des champs'
2001	do ig = 1, knon
2002	cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) &
2003	& + swdown(ig) / FLOAT(nexca)
2004	cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) &
2005	& + (lwdown(ig) + fluxlat(ig) +fluxsens(ig))&
2006	& / FLOAT(nexca)
2007	cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) &
2008	& + precip_rain(ig) / FLOAT(nexca)
2009	cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) &
2010	& + precip_snow(ig) / FLOAT(nexca)
2011	cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) &
2012	& + evap(ig) / FLOAT(nexca)
2013	cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) &
2014	& + tsurf(ig) / FLOAT(nexca)
2015	cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) &
2016	& + fder(ig) / FLOAT(nexca)
2017	cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) &
2018	& + albsol(ig) / FLOAT(nexca)
2019	cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) &
2020	& + taux(ig) / FLOAT(nexca)
2021	cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) &
2022	& + tauy(ig) / FLOAT(nexca)
2023	! -- LOOP
2024	IF (cpl_index .EQ. 1) THEN
2025	cpl_windsp(ig,cpl_index) = cpl_windsp(ig,cpl_index) &
2026	& + windsp(ig) / FLOAT(nexca)
2027	ENDIF
2028	! -- LOOP
2029	enddo
2030	IF (cpl_index .EQ. 1) THEN
2031	cpl_rriv(:,:) = cpl_rriv(:,:) + tmp_rriv(:,:) / FLOAT(nexca)
2032	cpl_rcoa(:,:) = cpl_rcoa(:,:) + tmp_rcoa(:,:) / FLOAT(nexca)
2033	cpl_rlic(:,:) = cpl_rlic(:,:) + tmp_rlic(:,:) / FLOAT(nexca)
2034	ENDIF
2035	endif
2036
2037	if (mod(itime, nexca) == 1) then
2038	!
2039	! Demande des champs au coupleur
2040	!
2041	! Si le domaine considere est l'ocean, on lit les champs venant du coupleur
2042	!
2043	if (nisurf == is_oce .and. .not. cumul) then
2044	if (check) write(,)'rentree fromcpl, itime-1 = ',itime-1
2045	#ifdef CPP_COUPLE
2046	#ifdef CPP_PSMILE
2047	il_time_secs=(itime-1)*dtime
2048	CALL fromcpl(il_time_secs, iim, jjphy_nb, &
2049	& read_sst, read_sic, read_sit, read_alb_sic)
2050	print *,read_sst
2051	#else
2052	if (.not. monocpu) then
2053	abort_message='coupleur parallele uniquement avec PSMILE'
2054	call abort_gcm(modname,abort_message,1)
2055	endif
2056
2057	call fromcpl(itime-1,(jjm+1)*iim, &
2058	& read_sst, read_sic, read_sit, read_alb_sic)
2059	#endif
2060	#endif
2061	!
2062	! sorties NETCDF des champs recus
2063	!
2064	!ym ndexcs(:)=0
2065	!ym itau_w = itau_phy + itime
2066	!ym CALL histwrite(nidcs,cl_read(1),itau_w,read_sst,iim*(jjm+1),ndexcs)
2067	!ym CALL histwrite(nidcs,cl_read(2),itau_w,read_sic,iim*(jjm+1),ndexcs)
2068	!ym CALL histwrite(nidcs,cl_read(3),itau_w,read_alb_sic,iim*(jjm+1),ndexcs)
2069	!ym CALL histwrite(nidcs,cl_read(4),itau_w,read_sit,iim*(jjm+1),ndexcs)
2070	!ym CALL histsync(nidcs)
2071	! pas utile IF (npas-itime.LT.nexca )CALL histclo(nidcs)
2072
2073	!ym do j = 1, jjm + 1
2074	do j = 1, jjphy_nb
2075	do ig = 1, iim
2076	if (abs(1. - read_sic(ig,j)) < 0.00001) then
2077	read_sst(ig,j) = RTT - 1.8
2078	read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
2079	read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
2080	else if (abs(read_sic(ig,j)) < 0.00001) then
2081	read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
2082	read_sit(ig,j) = read_sst(ig,j)
2083	read_alb_sic(ig,j) = 0.6
2084	else
2085	read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
2086	read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
2087	read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
2088	endif
2089	enddo
2090	enddo
2091	!
2092	! transformer read_sic en pctsrf_sav
2093	!
2094	call cpl2gath(read_sic, tamp_sic , klon, klon,iim,jjphy_nb, unity)
2095	do ig = 1, klon
2096	IF (pctsrf(ig,is_oce) > epsfra .OR. &
2097	& pctsrf(ig,is_sic) > epsfra) THEN
2098	pctsrf_sav(ig,is_sic) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
2099	& * tamp_sic(ig)
2100	pctsrf_sav(ig,is_oce) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
2101	& - pctsrf_sav(ig,is_sic)
2102	endif
2103	enddo
2104	!
2105	! Pour rattraper des erreurs d'arrondis
2106	!
2107	where (abs(pctsrf_sav(:,is_sic)) .le. 2.*epsilon(pctsrf_sav(1,is_sic)))
2108	pctsrf_sav(:,is_sic) = 0.
2109	pctsrf_sav(:,is_oce) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
2110	endwhere
2111	where (abs(pctsrf_sav(:,is_oce)) .le. 2.*epsilon(pctsrf_sav(1,is_oce)))
2112	pctsrf_sav(:,is_sic) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
2113	pctsrf_sav(:,is_oce) = 0.
2114	endwhere
2115	if (minval(pctsrf_sav(:,is_oce)) < 0.) then
2116	write(,)'Pb fraction ocean inferieure a 0'
2117	write(,)'au point ',minloc(pctsrf_sav(:,is_oce))
2118	write(,)'valeur = ',minval(pctsrf_sav(:,is_oce))
2119	abort_message = 'voir ci-dessus'
2120	call abort_gcm(modname,abort_message,1)
2121	endif
2122	if (minval(pctsrf_sav(:,is_sic)) < 0.) then
2123	write(,)'Pb fraction glace inferieure a 0'
2124	write(,)'au point ',minloc(pctsrf_sav(:,is_sic))
2125	write(,)'valeur = ',minval(pctsrf_sav(:,is_sic))
2126	abort_message = 'voir ci-dessus'
2127	call abort_gcm(modname,abort_message,1)
2128	endif
2129	endif
2130	endif ! fin mod(itime, nexca) == 1
2131
2132	if (mod(itime, nexca) == 0) then
2133	!
2134	! allocation memoire
2135	if (nisurf == is_oce .and. (.not. cumul) ) then
2136	sum_error = 0
2137	allocate(tmp_sols(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2138	allocate(tmp_nsol(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2139	allocate(tmp_rain(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2140	allocate(tmp_snow(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2141	allocate(tmp_evap(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2142	allocate(tmp_tsol(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2143	allocate(tmp_fder(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2144	allocate(tmp_albe(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2145	allocate(tmp_taux(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2146	allocate(tmp_tauy(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2147	! -- LOOP
2148	allocate(tmp_windsp(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2149	! -- LOOP
2150	!!$ allocate(tmp_rriv(iim,jjm+1,2), stat=error); sum_error = sum_error + error
2151	!!$ allocate(tmp_rcoa(iim,jjm+1,2), stat=error); sum_error = sum_error + error
2152	if (sum_error /= 0) then
2153	abort_message='Pb allocation variables couplees pour l''ecriture'
2154	call abort_gcm(modname,abort_message,1)
2155	endif
2156	endif
2157
2158	!
2159	! Mise sur la bonne grille des champs a passer au coupleur
2160	!
2161	cpl_index = 1
2162	if (nisurf == is_sic) cpl_index = 2
2163	call gath2cpl(cpl_sols(1,cpl_index), tmp_sols(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2164	call gath2cpl(cpl_nsol(1,cpl_index), tmp_nsol(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2165	call gath2cpl(cpl_rain(1,cpl_index), tmp_rain(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2166	call gath2cpl(cpl_snow(1,cpl_index), tmp_snow(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2167	call gath2cpl(cpl_evap(1,cpl_index), tmp_evap(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2168	call gath2cpl(cpl_tsol(1,cpl_index), tmp_tsol(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2169	call gath2cpl(cpl_fder(1,cpl_index), tmp_fder(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2170	call gath2cpl(cpl_albe(1,cpl_index), tmp_albe(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2171	call gath2cpl(cpl_taux(1,cpl_index), tmp_taux(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2172	call gath2cpl(cpl_tauy(1,cpl_index), tmp_tauy(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2173	! -- LOOP
2174	call gath2cpl(cpl_windsp(1,cpl_index), tmp_windsp(1,1,cpl_index), klon, knon,iim,jjphy_nb, knindex)
2175	! -- LOOP
2176
2177	!
2178	! Si le domaine considere est la banquise, on envoie les champs au coupleur
2179	!
2180	if (nisurf == is_sic .and. cumul) then
2181	wri_rain = 0.; wri_snow = 0.; wri_rcoa = 0.; wri_rriv = 0.
2182	wri_taux = 0.; wri_tauy = 0.
2183	! -- LOOP
2184	wri_windsp = 0.
2185	! -- LOOP
2186	call gath2cpl(pctsrf(1,is_oce), tamp_srf(1,1,1), klon, klon, iim, jjphy_nb, tamp_ind)
2187	call gath2cpl(pctsrf(1,is_sic), tamp_srf(1,1,2), klon, klon, iim, jjphy_nb, tamp_ind)
2188
2189	wri_sol_ice = tmp_sols(:,:,2)
2190	wri_sol_sea = tmp_sols(:,:,1)
2191	wri_nsol_ice = tmp_nsol(:,:,2)
2192	wri_nsol_sea = tmp_nsol(:,:,1)
2193	wri_fder_ice = tmp_fder(:,:,2)
2194	wri_evap_ice = tmp_evap(:,:,2)
2195	wri_evap_sea = tmp_evap(:,:,1)
2196	! -- LOOP
2197	wri_windsp = tmp_windsp(:,:,1)
2198	! -- LOOP
2199
2200	!!$PB
2201	wri_rriv = cpl_rriv(:,:)
2202	wri_rcoa = cpl_rcoa(:,:)
2203
2204	!ym !! ATTENTION ICI
2205
2206	!ym DO j = 1, jjm + 1
2207	!ym wri_calv(:,j) = sum(cpl_rlic(:,j)) / iim
2208	!ym enddo
2209
2210	!Essai OM+JLD : ca marche !!!! (17 mars 2006)
2211	tamp_srf(:,:,3)=0.
2212	CALL gath2cpl( pctsrf(1,is_lic), tamp_srf(1,1,3), klon, klon, iim, jjphy_nb, tamp_ind)
2213
2214	!YM pour retrouver resultat avant tamp_srf(:,3)=1.
2215
2216	DO j = 1, jjphy_nb
2217	wri_calv(:,j) = DOT_PRODUCT (cpl_rlic(1:iim,j), tamp_srf(1:iim,j,3)) / REAL(iim)
2218	ENDDO
2219
2220	!ym wri_calv(:,:)=0.
2221	!ym DO j = 1, jjphy_nb
2222	!ym wri_calv(:,j) = sum(cpl_rlic(:,j))/iim
2223	!ym enddo
2224
2225	IF (.NOT. monocpu) THEN
2226	if (phy_rank /= 0) then
2227	#ifdef CPP_PARA
2228	call MPI_RECV(Up,1,MPI_REAL8,phy_rank-1,1234,COMM_LMDZ,status,ierr)
2229	call MPI_SEND(wri_calv(1,1),1,MPI_REAL8,phy_rank-1,1234,COMM_LMDZ,ierr)
2230	#endif
2231	endif
2232
2233	if (phy_rank /= phy_size-1) then
2234	#ifdef CPP_PARA
2235	call MPI_SEND(wri_calv(1,jjphy_nb),1,MPI_REAL8,phy_rank+1,1234,COMM_LMDZ,ierr)
2236	call MPI_RECV(down,1,MPI_REAL8,phy_rank+1,1234,COMM_LMDZ,status,ierr)
2237	#endif
2238	endif
2239
2240	if (phy_rank /=0 .and. iiphy_begin /=1) then
2241	Up=Up+wri_calv(iim,1)
2242	wri_calv(:,1)=Up
2243	endif
2244
2245	if (phy_rank /=phy_size-1 .and. iiphy_end /= iim) then
2246	Down=Down+wri_calv(1,jjphy_nb)
2247	wri_calv(:,jjphy_nb)=Down
2248	endif
2249	ENDIF
2250
2251	where (tamp_zmasq /= 1.)
2252	deno = tamp_srf(:,:,1) + tamp_srf(:,:,2)
2253	wri_rain = tmp_rain(:,:,1) * tamp_srf(:,:,1) / deno + &
2254	& tmp_rain(:,:,2) * tamp_srf(:,:,2) / deno
2255	wri_snow = tmp_snow(:,:,1) * tamp_srf(:,:,1) / deno + &
2256	& tmp_snow(:,:,2) * tamp_srf(:,:,2) / deno
2257	wri_taux = tmp_taux(:,:,1) * tamp_srf(:,:,1) / deno + &
2258	& tmp_taux(:,:,2) * tamp_srf(:,:,2) / deno
2259	wri_tauy = tmp_tauy(:,:,1) * tamp_srf(:,:,1) / deno + &
2260	& tmp_tauy(:,:,2) * tamp_srf(:,:,2) / deno
2261	endwhere
2262	!
2263	! pour simuler la fonte des glaciers antarctiques
2264	!
2265	!$$$ wri_rain = wri_rain &
2266	!$$$ & + coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)
2267	! wri_calv = coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)
2268	!
2269	! on passe les coordonn�s de la grille
2270	!
2271
2272	!ym CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,tmp_lon)
2273	!ym CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,tmp_lat)
2274
2275	CALL phy2dyn(rlon,tmp_lon,1)
2276	CALL phy2dyn(rlat,tmp_lat,1)
2277
2278	!ym DO i = 1, iim
2279	!ym tmp_lon(i,1) = rlon(i+1)
2280	!ym tmp_lon(i,jjm + 1) = rlon(i+1)
2281	!ym ENDDO
2282
2283	!
2284	! sortie netcdf des champs pour le changement de repere
2285	!
2286	IF (monocpu) THEN
2287	ndexct(:)=0
2288	CALL histwrite(nidct,'tauxe',itau_w,wri_taux,iim*(jjm+1),ndexct)
2289	CALL histwrite(nidct,'tauyn',itau_w,wri_tauy,iim*(jjm+1),ndexct)
2290	CALL histwrite(nidct,'tmp_lon',itau_w,tmp_lon,iim*(jjm+1),ndexct)
2291	CALL histwrite(nidct,'tmp_lat',itau_w,tmp_lat,iim*(jjm+1),ndexct)
2292	ENDIF
2293	!
2294	! calcul 3 coordonn�s du vent
2295	!
2296	CALL atm2geo (iim , jjphy_nb, wri_taux, wri_tauy, tmp_lon, tmp_lat, &
2297	& wri_tauxx, wri_tauyy, wri_tauzz )
2298	!
2299	! sortie netcdf des champs apres changement de repere et juste avant
2300	! envoi au coupleur
2301	!
2302	IF (monocpu) THEN
2303	CALL histwrite(nidct,cl_writ(8),itau_w,wri_sol_ice,iim*(jjm+1),ndexct)
2304	CALL histwrite(nidct,cl_writ(9),itau_w,wri_sol_sea,iim*(jjm+1),ndexct)
2305	CALL histwrite(nidct,cl_writ(10),itau_w,wri_nsol_ice,iim*(jjm+1),ndexct)
2306	CALL histwrite(nidct,cl_writ(11),itau_w,wri_nsol_sea,iim*(jjm+1),ndexct)
2307	CALL histwrite(nidct,cl_writ(12),itau_w,wri_fder_ice,iim*(jjm+1),ndexct)
2308	CALL histwrite(nidct,cl_writ(13),itau_w,wri_evap_ice,iim*(jjm+1),ndexct)
2309	CALL histwrite(nidct,cl_writ(14),itau_w,wri_evap_sea,iim*(jjm+1),ndexct)
2310	CALL histwrite(nidct,cl_writ(15),itau_w,wri_rain,iim*(jjm+1),ndexct)
2311	CALL histwrite(nidct,cl_writ(16),itau_w,wri_snow,iim*(jjm+1),ndexct)
2312	CALL histwrite(nidct,cl_writ(17),itau_w,wri_rcoa,iim*(jjm+1),ndexct)
2313	CALL histwrite(nidct,cl_writ(18),itau_w,wri_rriv,iim*(jjm+1),ndexct)
2314	CALL histwrite(nidct,cl_writ(19),itau_w,wri_calv,iim*(jjm+1),ndexct)
2315	CALL histwrite(nidct,cl_writ(1),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
2316	CALL histwrite(nidct,cl_writ(2),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
2317	CALL histwrite(nidct,cl_writ(3),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
2318	CALL histwrite(nidct,cl_writ(4),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
2319	CALL histwrite(nidct,cl_writ(5),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
2320	CALL histwrite(nidct,cl_writ(6),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
2321	! -- LOOP
2322	CALL histwrite(nidct,cl_writ(7),itau_w,wri_windsp,iim*(jjm+1),ndexct)
2323	! -- LOOP
2324	CALL histsync(nidct)
2325	ENDIF
2326	! pas utile IF (lafin) CALL histclo(nidct)
2327	#ifdef CPP_COUPLE
2328	#ifdef CPP_PSMILE
2329	il_time_secs=(itime-1)*dtime
2330
2331	CALL intocpl(il_time_secs, iim, jjm+1, wri_sol_ice, wri_sol_sea, wri_nsol_ice,&
2332	& wri_nsol_sea, wri_fder_ice, wri_evap_ice, wri_evap_sea, wri_rain, &
2333	& wri_snow, wri_rcoa, wri_rriv, wri_calv, wri_tauxx, wri_tauyy, &
2334	& wri_tauzz, wri_tauxx, wri_tauyy, wri_tauzz, &
2335	! -- LOOP
2336	& wri_windsp,lafin)
2337	! -- LOOP
2338	#else
2339	call intocpl(itime, (jjm+1)*iim, wri_sol_ice, wri_sol_sea, wri_nsol_ice,&
2340	& wri_nsol_sea, wri_fder_ice, wri_evap_ice, wri_evap_sea, wri_rain, &
2341	& wri_snow, wri_rcoa, wri_rriv, wri_calv, wri_tauxx, wri_tauyy, &
2342	& wri_tauzz, wri_tauxx, wri_tauyy, wri_tauzz, &
2343	! -- LOOP
2344	& wri_windsp,lafin)
2345	! -- LOOP
2346	#endif
2347	#endif
2348	!
2349	cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
2350	cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
2351	cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.; cpl_rlic = 0.
2352	! -- LOOP
2353	cpl_windsp = 0.
2354	! -- LOOP
2355	!
2356	! deallocation memoire variables temporaires
2357	!
2358	sum_error = 0
2359	deallocate(tmp_sols, stat=error); sum_error = sum_error + error
2360	deallocate(tmp_nsol, stat=error); sum_error = sum_error + error
2361	deallocate(tmp_rain, stat=error); sum_error = sum_error + error
2362	deallocate(tmp_snow, stat=error); sum_error = sum_error + error
2363	deallocate(tmp_evap, stat=error); sum_error = sum_error + error
2364	deallocate(tmp_fder, stat=error); sum_error = sum_error + error
2365	deallocate(tmp_tsol, stat=error); sum_error = sum_error + error
2366	deallocate(tmp_albe, stat=error); sum_error = sum_error + error
2367	deallocate(tmp_taux, stat=error); sum_error = sum_error + error
2368	deallocate(tmp_tauy, stat=error); sum_error = sum_error + error
2369	! -- LOOP
2370	deallocate(tmp_windsp, stat=error); sum_error = sum_error + error
2371	! -- LOOP
2372	!!$PB
2373	!!$ deallocate(tmp_rriv, stat=error); sum_error = sum_error + error
2374	!!$ deallocate(tmp_rcoa, stat=error); sum_error = sum_error + error
2375	if (sum_error /= 0) then
2376	abort_message='Pb deallocation variables couplees'
2377	call abort_gcm(modname,abort_message,1)
2378	endif
2379
2380	endif
2381
2382	endif ! fin (mod(itime, nexca) == 0)
2383	!
2384	! on range les variables lues/sauvegardees dans les bonnes variables de sortie
2385	!
2386	if (nisurf == is_oce) then
2387	call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjphy_nb, knindex)
2388	else if (nisurf == is_sic) then
2389	call cpl2gath(read_sit, tsurf_new, klon, knon,iim,jjphy_nb, knindex)
2390	call cpl2gath(read_alb_sic, alb_new, klon, knon,iim,jjphy_nb, knindex)
2391	endif
2392	pctsrf_new(:,nisurf) = pctsrf_sav(:,nisurf)
2393
2394	if (mod(itime, nexca) == -1) then
2395	tmp_field=0.
2396	do i = 1, knon
2397	ig = knindex(i)
2398	tmp_field(ig) = 1.
2399	enddo
2400	call WriteField_phy('knindex',tmp_field,1)
2401
2402	tmp_field=0.
2403	do i = 1, knon
2404	ig = knindex(i)
2405	tmp_field(ig) = tsurf_new(i)
2406	enddo
2407	call WriteField_phy('tsurf_new',tmp_field,1)
2408
2409	tmp_field=0.
2410	do i = 1, knon
2411	ig = knindex(i)
2412	tmp_field(ig) = alb_new(i)
2413	enddo
2414	call WriteField_phy('alb_new',tmp_field,1)
2415
2416	! tmp_field=0.
2417	! do i = 1, knon
2418	! ig = knindex(i)
2419	! tmp_field(ig) = pctsrf_new(i,nisurf)
2420	! enddo
2421	call WriteField_phy('pctsrf_new', pctsrf_new(:,nisurf),1)
2422	endif
2423	!ym do j=1,jjphy_nb
2424	!ym do i=1,iim
2425	!ym print *,phy_rank,'read_sst(',i,',',j,')=',read_sst(i,j)
2426	!ym enddo
2427	!ym enddo
2428
2429	!ym do i=1,knon
2430	!ym print *,phy_rank,'tsurf_new(',i,')=',tsurf_new(i)
2431	!ym enddo
2432	! if (lafin) call quitcpl
2433
2434	END SUBROUTINE interfoce_cpl
2435	!
2436	!#########################################################################
2437	!
2438	SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &
2439	& radsol, fluxo, fluxg, pctsrf, &
2440	& tslab, seaice, pctsrf_slab)
2441	!
2442	! Cette routine calcule la temperature d'un slab ocean, la glace de mer
2443	! et les pourcentages de la maille couverte par l'ocean libre et/ou
2444	! la glace de mer pour un "slab" ocean de 50m
2445	!
2446	! Conception: Laurent Li
2447	! Re-ecriture + adaptation LMDZ4: I. Musat
2448	!
2449	! input:
2450	! klon nombre total de points de grille
2451	! debut logical: 1er appel a la physique
2452	! itap numero du pas de temps
2453	! dtime pas de temps de la physique (en s)
2454	! ijour jour dans l'annee en cours
2455	! radsol rayonnement net au sol (LW + SW)
2456	! fluxo flux turbulent (sensible + latent) sur les mailles oceaniques
2457	! fluxg flux de conduction entre la surface de la glace de mer et l'ocean
2458	! pctsrf tableau des pourcentages de surface de chaque maille
2459	! output:
2460	! tslab temperature de l'ocean libre
2461	! seaice glace de mer (kg/m2)
2462	! pctsrf_slab "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
2463	!
2464	#include "indicesol.inc"
2465	#include "clesphys.inc"
2466	! Parametres d'entree
2467	integer, intent(IN) :: klon
2468	logical, intent(IN) :: debut
2469	INTEGER, intent(IN) :: itap
2470	REAL, intent(IN) :: dtime
2471	INTEGER, intent(IN) :: ijour
2472	REAL, dimension(klon), intent(IN) :: radsol
2473	REAL, dimension(klon), intent(IN) :: fluxo
2474	REAL, dimension(klon), intent(IN) :: fluxg
2475	real, dimension(klon, nbsrf), intent(IN) :: pctsrf
2476	! Parametres de sortie
2477	real, dimension(klon), intent(INOUT) :: tslab
2478	real, dimension(klon), intent(INOUT) :: seaice ! glace de mer (kg/m2)
2479	real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_slab
2480	!
2481	! Variables locales :
2482	REAL :: amn, amx
2483	INTEGER, save :: lmt_pas, julien, idayvrai
2484	REAL, parameter :: unjour=86400.
2485	real, allocatable, dimension(:), save :: tmp_tslab, tmp_seaice
2486	REAL, allocatable, dimension(:), save :: slab_bils
2487	REAL, allocatable, dimension(:), save :: lmt_bils
2488	logical,save :: check = .false.
2489	!
2490	REAL, parameter :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
2491	REAL, parameter :: cbing=0.334e+05 ! J/kg
2492	real, dimension(klon) :: siceh !hauteur de la glace de mer (m)
2493	INTEGER :: i
2494	integer :: sum_error, error
2495	REAL :: zz, za, zb
2496	!
2497	character (len = 80) :: abort_message
2498	character (len = 20) :: modname = 'interfoce_slab'
2499	!
2500	julien = MOD(ijour,360)
2501	sum_error = 0
2502	IF (debut) THEN
2503	allocate(slab_bils(klon), stat = error); sum_error = sum_error + error
2504	allocate(lmt_bils(klon), stat = error); sum_error = sum_error + error
2505	allocate(tmp_tslab(klon), stat = error); sum_error = sum_error + error
2506	allocate(tmp_seaice(klon), stat = error); sum_error = sum_error + error
2507	if (sum_error /= 0) then
2508	abort_message='Pb allocation var. slab_bils,lmt_bils,tmp_tslab,tmp_seaice'
2509	call abort_gcm(modname,abort_message,1)
2510	endif
2511	tmp_tslab=tslab
2512	tmp_seaice=seaice
2513	lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
2514	!
2515	IF (check) THEN
2516	PRINT*,'interfoce_slab klon, debut, itap, dtime, ijour, &
2517	& lmt_pas ', klon, debut, itap, dtime, ijour, &
2518	& lmt_pas
2519	ENDIF !check
2520	!
2521	PRINT, '***********************'
2522	PRINT*, 'SLAB OCEAN est actif, prenez precautions !'
2523	PRINT, '***********************'
2524	!
2525	! a mettre un slab_bils aussi en force !!!
2526	!
2527	DO i = 1, klon
2528	slab_bils(i) = 0.0
2529	ENDDO
2530	!
2531	ENDIF !debut
2532	!
2533	IF (check ) THEN
2534	amn=MIN(tmp_tslab(1),1000.)
2535	amx=MAX(tmp_tslab(1),-1000.)
2536	DO i=2, klon
2537	amn=MIN(tmp_tslab(i),amn)
2538	amx=MAX(tmp_tslab(i),amx)
2539	ENDDO
2540	!
2541	PRINT*,' debut min max tslab',amn,amx
2542	!
2543	!!
2544	PRINT*,' itap,lmt_pas unjour',itap,lmt_pas,unjour
2545	ENDIF
2546	!!
2547	!
2548	pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)
2549	!
2550	! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
2551	!
2552	IF (MOD(itap,lmt_pas) .EQ. 1) THEN !1er pas de temps de la journee
2553	idayvrai = ijour
2554	CALL condsurf(julien,idayvrai, lmt_bils)
2555	ENDIF !(MOD(itap-1,lmt_pas) .EQ. 0) THEN
2556
2557	DO i = 1, klon
2558	IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
2559	& (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
2560	!
2561	! fabriquer de la glace si congelation atteinte:
2562	!
2563	IF (tmp_tslab(i).LT.(RTT-1.8)) THEN
2564	zz = (RTT-1.8)-tmp_tslab(i)
2565	tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz
2566	seaice(i) = tmp_seaice(i)
2567	tmp_tslab(i) = RTT-1.8
2568	ENDIF
2569	!
2570	! faire fondre de la glace si temperature est superieure a 0:
2571	!
2572	IF ((tmp_tslab(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN
2573	zz = cyang/cbing * (tmp_tslab(i)-RTT)
2574	zz = MIN(zz,tmp_seaice(i))
2575	tmp_seaice(i) = tmp_seaice(i) - zz
2576	seaice(i) = tmp_seaice(i)
2577	tmp_tslab(i) = tmp_tslab(i) - zz*cbing/cyang
2578	ENDIF
2579	!
2580	! limiter la glace de mer a 10 metres (10000 kg/m2)
2581	!
2582	IF(tmp_seaice(i).GT.45.) THEN
2583	tmp_seaice(i) = MIN(tmp_seaice(i),10000.0)
2584	ELSE
2585	tmp_seaice(i) = 0.
2586	ENDIF
2587	seaice(i) = tmp_seaice(i)
2588	siceh(i)=tmp_seaice(i)/1000. !en metres
2589	!
2590	! determiner la nature du sol (glace de mer ou ocean libre):
2591	!
2592	! on fait dependre la fraction de seaice "pctsrf(i,is_sic)"
2593	! de l'epaisseur de seaice :
2594	! pctsrf(i,is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm
2595	! et pctsrf(i,is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur
2596	!
2597
2598	IF(.NOT.ok_slab_sicOBS) then
2599	pctsrf_slab(i,is_sic)=MIN(siceh(i)/0.20, &
2600	& 1.-(pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)))
2601	pctsrf_slab(i,is_oce)=1.0 - &
2602	& (pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)+pctsrf_slab(i,is_sic))
2603	ELSE
2604	IF (i.EQ.1) print*,'cas ok_slab_sicOBS TRUE : passe sur la modif.'
2605	ENDIF !(.NOT.ok_slab_sicOBS) then
2606	ENDIF !pctsrf
2607	ENDDO
2608	!
2609	! Calculer le bilan du flux de chaleur au sol :
2610	!
2611	DO i = 1, klon
2612	za = radsol(i) + fluxo(i)
2613	zb = fluxg(i)
2614	IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
2615	& (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
2616	slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i,is_oce) &
2617	& +zb*pctsrf_slab(i,is_sic))/ FLOAT(lmt_pas)
2618	ENDIF
2619	ENDDO !klon
2620	!
2621	! calcul tslab
2622	!
2623	IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee
2624	!
2625	! calcul tslab
2626	amn=MIN(tmp_tslab(1),1000.)
2627	amx=MAX(tmp_tslab(1),-1000.)
2628	DO i = 1, klon
2629	IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
2630	& (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
2631	tmp_tslab(i) = tmp_tslab(i) + &
2632	& (slab_bils(i)-lmt_bils(i)) &
2633	& /cyang*unjour
2634	! on remet l'accumulation a 0
2635	slab_bils(i) = 0.
2636	ENDIF !pctsrf
2637	!
2638	IF (check) THEN
2639	IF(i.EQ.1) THEN
2640	PRINT*,'i tmp_tslab MOD(itap,lmt_pas).EQ.0 sicINTER',i,itap, &
2641	& tmp_tslab(i), tmp_seaice(i)
2642	ENDIF
2643	!
2644	amn=MIN(tmp_tslab(i),amn)
2645	amx=MAX(tmp_tslab(i),amx)
2646	ENDIF
2647	ENDDO !klon
2648	ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN
2649	!
2650	IF ( check ) THEN
2651	PRINT*,'fin min max tslab',amn,amx
2652	ENDIF
2653	!
2654	tslab = tmp_tslab
2655	seaice = tmp_seaice
2656	END SUBROUTINE interfoce_slab
2657	!
2658	!#########################################################################
2659	!
2660	SUBROUTINE interfoce_lim(itime, dtime, jour, &
2661	& klon_xx, nisurf, knon, knindex, &
2662	& debut, &
2663	& lmt_sst_p, pctsrf_new_p)
2664
2665	USE dimphy,klon=>klon2,klon2=>klon
2666
2667	#include "indicesol.inc"
2668
2669	! Cette routine sert d'interface entre le modele atmospherique et un fichier
2670	! de conditions aux limites
2671	!
2672	! L. Fairhead 02/2000
2673	!
2674	! input:
2675	! itime numero du pas de temps courant
2676	! dtime pas de temps de la physique (en s)
2677	! jour jour a lire dans l'annee
2678	! nisurf index de la surface a traiter (1 = sol continental)
2679	! knon nombre de points dans le domaine a traiter
2680	! knindex index des points de la surface a traiter
2681	! klon taille de la grille
2682	! debut logical: 1er appel a la physique (initialisation)
2683	!
2684	! output:
2685	! lmt_sst SST lues dans le fichier de CL
2686	! pctsrf_new sous-maille fractionnelle
2687	!
2688
2689
2690	! Parametres d'entree
2691	integer, intent(IN) :: itime
2692	real , intent(IN) :: dtime
2693	integer, intent(IN) :: jour
2694	integer, intent(in) :: klon_xx
2695	integer, intent(IN) :: nisurf
2696	integer, intent(IN) :: knon
2697	integer, dimension(klon2), intent(in) :: knindex
2698	logical, intent(IN) :: debut
2699
2700	! Parametres de sortie
2701	real, intent(out), dimension(klon2) :: lmt_sst_p
2702	real, intent(out), dimension(klon2,nbsrf) :: pctsrf_new_p
2703
2704	! real, dimension(klon) :: lmt_sst
2705	real, dimension(klon,nbsrf) :: pctsrf_new
2706
2707	! Variables locales
2708	integer :: ii
2709	INTEGER,save :: lmt_pas ! frequence de lecture des conditions limites
2710	! (en pas de physique)
2711	!$OMP THREADPRIVATE(lmt_pas)
2712	logical,save :: deja_lu ! pour indiquer que le jour a lire a deja
2713	! lu pour une surface precedente
2714	!$OMP THREADPRIVATE(deja_lu)
2715	integer,save :: jour_lu
2716	!$OMP THREADPRIVATE(jour_lu)
2717	integer :: ierr
2718	character (len = 20) :: modname = 'interfoce_lim'
2719	character (len = 80) :: abort_message
2720	character (len = 20),save :: fich ='limit.nc'
2721	!$OMP THREADPRIVATE(fich)
2722	logical, save :: newlmt = .TRUE.
2723	!$OMP THREADPRIVATE(newlmt)
2724	logical, save :: check = .FALSE.
2725	!$OMP THREADPRIVATE(check)
2726	! Champs lus dans le fichier de CL
2727	real, allocatable , save, dimension(:) :: sst_lu_p
2728	!$OMP THREADPRIVATE(sst_lu_p)
2729	real, allocatable , save, dimension(:) :: sst_lu_mpi
2730
2731	real, allocatable , save, dimension(:,:) :: pct_tmp_p
2732	!$OMP THREADPRIVATE(pct_tmp_p)
2733	real, allocatable , save, dimension(:,:) :: pct_tmp_mpi
2734	real, dimension(klon,nbsrf) :: pct_tmp
2735	real, dimension(klon) :: sst_lu
2736	real, dimension(klon) :: nat_lu
2737	!
2738	! quelques variables pour netcdf
2739	!
2740	#include "netcdf.inc"
2741	integer :: nid, nvarid
2742	integer, dimension(2) :: start, epais
2743	!
2744	! Fin d�laration
2745	!
2746
2747	if (debut .and. .not. allocated(sst_lu_p)) then
2748	lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
2749	jour_lu = jour - 1
2750	allocate(sst_lu_p(klon_omp))
2751	allocate(pct_tmp_p(klon_omp,nbsrf))
2752	endif
2753
2754	if ((jour - jour_lu) /= 0) deja_lu = .false.
2755
2756	if (check) write(,)modname,' :: jour, jour_lu, deja_lu', jour, jour_lu, deja_lu
2757	if (check) write(,)modname,' :: itime, lmt_pas ', itime, lmt_pas,dtime
2758
2759	! Tester d'abord si c'est le moment de lire le fichier
2760	if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then
2761	!
2762	! Ouverture du fichier
2763	!
2764	!$OMP MASTER
2765	if (.not. allocated(sst_lu_mpi)) allocate(sst_lu_mpi(klon_mpi))
2766	if (.not. allocated(pct_tmp_mpi)) allocate(pct_tmp_mpi(klon_mpi,nbsrf))
2767
2768	if (phy_rank==0) then
2769
2770	fich = trim(fich)
2771	ierr = NF_OPEN (fich, NF_NOWRITE,nid)
2772	if (ierr.NE.NF_NOERR) then
2773	abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
2774	call abort_gcm(modname,abort_message,1)
2775	endif
2776	!
2777	! La tranche de donnees a lire:
2778	!
2779	start(1) = 1
2780	start(2) = jour
2781	epais(1) = klon
2782	epais(2) = 1
2783	!
2784	if (newlmt) then
2785	!
2786	! Fraction "ocean"
2787	!
2788	ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)
2789	if (ierr /= NF_NOERR) then
2790	abort_message = 'Le champ <FOCE> est absent'
2791	call abort_gcm(modname,abort_message,1)
2792	endif
2793	#ifdef NC_DOUBLE
2794	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_oce))
2795	#else
2796	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_oce))
2797	#endif
2798	if (ierr /= NF_NOERR) then
2799	abort_message = 'Lecture echouee pour <FOCE>'
2800	call abort_gcm(modname,abort_message,1)
2801	endif
2802	!
2803	! Fraction "glace de mer"
2804	!
2805	ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)
2806	if (ierr /= NF_NOERR) then
2807	abort_message = 'Le champ <FSIC> est absent'
2808	call abort_gcm(modname,abort_message,1)
2809	endif
2810	#ifdef NC_DOUBLE
2811	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_sic))
2812	#else
2813	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_sic))
2814	#endif
2815	if (ierr /= NF_NOERR) then
2816	abort_message = 'Lecture echouee pour <FSIC>'
2817	call abort_gcm(modname,abort_message,1)
2818	endif
2819	!
2820	! Fraction "terre"
2821	!
2822	ierr = NF_INQ_VARID(nid, 'FTER', nvarid)
2823	if (ierr /= NF_NOERR) then
2824	abort_message = 'Le champ <FTER> est absent'
2825	call abort_gcm(modname,abort_message,1)
2826	endif
2827	#ifdef NC_DOUBLE
2828	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_ter))
2829	#else
2830	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_ter))
2831	#endif
2832	if (ierr /= NF_NOERR) then
2833	abort_message = 'Lecture echouee pour <FTER>'
2834	call abort_gcm(modname,abort_message,1)
2835	endif
2836	!
2837	! Fraction "glacier terre"
2838	!
2839	ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)
2840	if (ierr /= NF_NOERR) then
2841	abort_message = 'Le champ <FLIC> est absent'
2842	call abort_gcm(modname,abort_message,1)
2843	endif
2844	#ifdef NC_DOUBLE
2845	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_lic))
2846	#else
2847	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_lic))
2848	#endif
2849	if (ierr /= NF_NOERR) then
2850	abort_message = 'Lecture echouee pour <FLIC>'
2851	call abort_gcm(modname,abort_message,1)
2852	endif
2853	!
2854	else ! on en est toujours a rnatur
2855	!
2856	ierr = NF_INQ_VARID(nid, 'NAT', nvarid)
2857	if (ierr /= NF_NOERR) then
2858	abort_message = 'Le champ <NAT> est absent'
2859	call abort_gcm(modname,abort_message,1)
2860	endif
2861	#ifdef NC_DOUBLE
2862	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, nat_lu)
2863	#else
2864	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, nat_lu)
2865	#endif
2866	if (ierr /= NF_NOERR) then
2867	abort_message = 'Lecture echouee pour <NAT>'
2868	call abort_gcm(modname,abort_message,1)
2869	endif
2870	!
2871	! Remplissage des fractions de surface
2872	! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice
2873	!
2874	pct_tmp = 0.0
2875	do ii = 1, klon
2876	pct_tmp(ii,nint(nat_lu(ii)) + 1) = 1.
2877	enddo
2878
2879	!
2880	! On se retrouve avec ocean en 1 et terre en 2 alors qu'on veut le contraire
2881	!
2882	pctsrf_new = pct_tmp
2883	pctsrf_new (:,2)= pct_tmp (:,1)
2884	pctsrf_new (:,1)= pct_tmp (:,2)
2885	pct_tmp = pctsrf_new
2886	endif ! fin test sur newlmt
2887	!
2888	! Lecture SST
2889	!
2890	ierr = NF_INQ_VARID(nid, 'SST', nvarid)
2891	if (ierr /= NF_NOERR) then
2892	abort_message = 'Le champ <SST> est absent'
2893	call abort_gcm(modname,abort_message,1)
2894	endif
2895	#ifdef NC_DOUBLE
2896	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, sst_lu)
2897	#else
2898	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, sst_lu)
2899	#endif
2900	if (ierr /= NF_NOERR) then
2901	abort_message = 'Lecture echouee pour <SST>'
2902	call abort_gcm(modname,abort_message,1)
2903	endif
2904
2905	!
2906	! Fin de lecture
2907	!
2908	ierr = NF_CLOSE(nid)
2909	endif ! phyrank
2910	!
2911	! Recopie des variables dans les champs de sortie
2912	!
2913	call ScatterField(sst_lu,sst_lu_mpi,1)
2914	call ScatterField(pct_tmp(:,is_oce),pct_tmp_mpi(:,is_oce),1)
2915	call ScatterField(pct_tmp(:,is_sic),pct_tmp_mpi(:,is_sic),1)
2916	!$OMP END MASTER
2917	!$OMP BARRIER
2918	call ScatterField_omp(sst_lu_mpi,sst_lu_p,1)
2919	call ScatterField_omp(pct_tmp_mpi(:,is_oce),pct_tmp_p(:,is_oce),1)
2920	call ScatterField_omp(pct_tmp_mpi(:,is_sic),pct_tmp_p(:,is_sic),1)
2921	deja_lu = .true.
2922	jour_lu = jour
2923	endif
2924
2925	lmt_sst_p = 999999999.
2926
2927	do ii = 1, knon
2928	lmt_sst_p(ii) = sst_lu_p(knindex(ii))
2929	enddo
2930
2931	do ii=1,klon2
2932	pctsrf_new_p(ii,is_oce)=pct_tmp_p(ii,is_oce)
2933	pctsrf_new_p(ii,is_sic)=pct_tmp_p(ii,is_sic)
2934	enddo
2935
2936
2937	END SUBROUTINE interfoce_lim
2938
2939	!
2940	!#########################################################################
2941	!
2942	SUBROUTINE interfsur_lim(itime, dtime, jour, &
2943	& klon_xx, nisurf, knon, knindex, &
2944	& debut, &
2945	& lmt_alb_p, lmt_rug_p)
2946
2947	USE dimphy,klon=>klon2,klon2=>klon
2948
2949	! Cette routine sert d'interface entre le modele atmospherique et un fichier
2950	! de conditions aux limites
2951	!
2952	! L. Fairhead 02/2000
2953	!
2954	! input:
2955	! itime numero du pas de temps courant
2956	! dtime pas de temps de la physique (en s)
2957	! jour jour a lire dans l'annee
2958	! nisurf index de la surface a traiter (1 = sol continental)
2959	! knon nombre de points dans le domaine a traiter
2960	! knindex index des points de la surface a traiter
2961	! klon taille de la grille
2962	! debut logical: 1er appel a la physique (initialisation)
2963	!
2964	! output:
2965	! lmt_sst SST lues dans le fichier de CL
2966	! lmt_alb Albedo lu
2967	! lmt_rug longueur de rugosit�lue
2968	! pctsrf_new sous-maille fractionnelle
2969	!
2970
2971
2972	! Parametres d'entree
2973	integer, intent(IN) :: itime
2974	real , intent(IN) :: dtime
2975	integer, intent(IN) :: jour
2976	integer, intent(IN) :: nisurf
2977	integer, intent(IN) :: knon
2978	integer, intent(IN) :: klon_xx
2979	integer, dimension(klon2), intent(in) :: knindex
2980	logical, intent(IN) :: debut
2981
2982	! Parametres de sortie
2983	real, intent(out), dimension(klon2) :: lmt_alb_p
2984	real, intent(out), dimension(klon2) :: lmt_rug_p
2985
2986	! real, dimension(klon) :: lmt_alb
2987	! real, dimension(klon) :: lmt_rug
2988
2989	! Variables locales
2990	integer :: ii
2991	integer,save :: lmt_pas ! frequence de lecture des conditions limites
2992	! (en pas de physique)
2993	!$OMP THREADPRIVATE(lmt_pas)
2994	logical,save :: deja_lu_sur! pour indiquer que le jour a lire a deja
2995	! lu pour une surface precedente
2996	!$OMP THREADPRIVATE(deja_lu_sur)
2997	integer,save :: jour_lu_sur
2998	!$OMP THREADPRIVATE(jour_lu_sur)
2999	integer :: ierr
3000	character (len = 20) :: modname = 'interfsur_lim'
3001	character (len = 80) :: abort_message
3002	character (len = 20),save :: fich ='limit.nc'
3003	!$OMP THREADPRIVATE(fich)
3004	logical,save :: newlmt = .false.
3005	!$OMP THREADPRIVATE(newlmt)
3006	logical,save :: check = .false.
3007	!$OMP THREADPRIVATE(check)
3008	! Champs lus dans le fichier de CL
3009	real, allocatable , save, dimension(:) :: alb_lu_p, rug_lu_p
3010	!$OMP THREADPRIVATE(alb_lu_p, rug_lu_p)
3011	real, allocatable , save, dimension(:) :: alb_lu_mpi, rug_lu_mpi
3012	real, dimension(klon) :: alb_lu, rug_lu
3013	!
3014	! quelques variables pour netcdf
3015	!
3016	#include "netcdf.inc"
3017	integer ,save :: nid, nvarid
3018	!$OMP THREADPRIVATE(nid, nvarid)
3019	integer, dimension(2),save :: start, epais
3020	!$OMP THREADPRIVATE(start, epais)
3021	!
3022	! Fin d�laration
3023	!
3024
3025	if (debut) then
3026	lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
3027	jour_lu_sur = jour - 1
3028	allocate(alb_lu_p(klon_omp))
3029	allocate(rug_lu_p(klon_omp))
3030	endif
3031
3032	if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false.
3033
3034	if (check) write(,)modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, deja_lu_sur
3035	if (check) write(,)modname,':: itime, lmt_pas', itime, lmt_pas
3036	if (check) call flush(6)
3037
3038	! Tester d'abord si c'est le moment de lire le fichier
3039	if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu_sur) then
3040
3041	!$OMP MASTER
3042	if (.not. allocated(alb_lu_mpi)) allocate(alb_lu_mpi(klon_mpi))
3043	if (.not. allocated(rug_lu_mpi)) allocate(rug_lu_mpi(klon_mpi))
3044	if (phy_rank==0) then
3045	!
3046	! Ouverture du fichier
3047	!
3048	fich = trim(fich)
3049	IF (check) WRITE(,)modname,' ouverture fichier ',fich
3050	if (check) CALL flush(6)
3051	ierr = NF_OPEN (fich, NF_NOWRITE,nid)
3052	if (ierr.NE.NF_NOERR) then
3053	abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
3054	call abort_gcm(modname,abort_message,1)
3055	endif
3056	!
3057	! La tranche de donnees a lire:
3058
3059	start(1) = 1
3060	start(2) = jour
3061	epais(1) = klon
3062	epais(2) = 1
3063	!
3064	! Lecture Albedo
3065	!
3066	ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
3067	if (ierr /= NF_NOERR) then
3068	abort_message = 'Le champ <ALB> est absent'
3069	call abort_gcm(modname,abort_message,1)
3070	endif
3071	#ifdef NC_DOUBLE
3072	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, alb_lu)
3073	#else
3074	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)
3075	#endif
3076	if (ierr /= NF_NOERR) then
3077	abort_message = 'Lecture echouee pour <ALB>'
3078	call abort_gcm(modname,abort_message,1)
3079	endif
3080	!
3081	! Lecture rugosit�!
3082	ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
3083	if (ierr /= NF_NOERR) then
3084	abort_message = 'Le champ <RUG> est absent'
3085	call abort_gcm(modname,abort_message,1)
3086	endif
3087	#ifdef NC_DOUBLE
3088	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, rug_lu)
3089	#else
3090	ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)
3091	#endif
3092	if (ierr /= NF_NOERR) then
3093	abort_message = 'Lecture echouee pour <RUG>'
3094	call abort_gcm(modname,abort_message,1)
3095	endif
3096
3097	!
3098	! Fin de lecture
3099	!
3100	ierr = NF_CLOSE(nid)
3101
3102
3103	endif !! phyrank
3104
3105	call ScatterField(alb_lu,alb_lu_mpi,1)
3106	call ScatterField(rug_lu,rug_lu_mpi,1)
3107	!$OMP END MASTER
3108	!$OMP BARRIER
3109
3110	call ScatterField_omp(alb_lu_mpi,alb_lu_p,1)
3111	call ScatterField_omp(rug_lu_mpi,rug_lu_p,1)
3112
3113	deja_lu_sur = .true.
3114	jour_lu_sur = jour
3115
3116
3117	endif
3118
3119	!
3120	! Recopie des variables dans les champs de sortie
3121	!
3122	!!$ lmt_alb(:) = 0.0
3123	!!$ lmt_rug(:) = 0.0
3124
3125	lmt_alb_p(:) = 999999.
3126	lmt_rug_p(:) = 999999.
3127	DO ii = 1, knon
3128	lmt_alb_p(ii) = alb_lu_p(knindex(ii))
3129	lmt_rug_p(ii) = rug_lu_p(knindex(ii))
3130	enddo
3131
3132
3133	END SUBROUTINE interfsur_lim
3134
3135	!
3136	!#########################################################################
3137	!
3138
3139	SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime, &
3140	& tsurf, p1lay, cal, beta, coef1lay, ps, &
3141	& precip_rain, precip_snow, snow, qsurf, &
3142	& radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
3143	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
3144	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
3145	USE dimphy,only : omp_rank
3146	! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
3147	! une temperature de surface (au cas ou ok_veget = false)
3148	!
3149	! L. Fairhead 4/2000
3150	!
3151	! input:
3152	! knon nombre de points a traiter
3153	! nisurf surface a traiter
3154	! tsurf temperature de surface
3155	! p1lay pression 1er niveau (milieu de couche)
3156	! cal capacite calorifique du sol
3157	! beta evap reelle
3158	! coef1lay coefficient d'echange
3159	! ps pression au sol
3160	! precip_rain precipitations liquides
3161	! precip_snow precipitations solides
3162	! snow champs hauteur de neige
3163	! runoff runoff en cas de trop plein
3164	! petAcoef coeff. A de la resolution de la CL pour t
3165	! peqAcoef coeff. A de la resolution de la CL pour q
3166	! petBcoef coeff. B de la resolution de la CL pour t
3167	! peqBcoef coeff. B de la resolution de la CL pour q
3168	! radsol rayonnement net aus sol (LW + SW)
3169	! dif_grnd coeff. diffusion vers le sol profond
3170	!
3171	! output:
3172	! tsurf_new temperature au sol
3173	! qsurf humidite de l'air au dessus du sol
3174	! fluxsens flux de chaleur sensible
3175	! fluxlat flux de chaleur latente
3176	! dflux_s derivee du flux de chaleur sensible / Ts
3177	! dflux_l derivee du flux de chaleur latente / Ts
3178	!
3179
3180	#include "YOETHF.inc"
3181	#include "FCTTRE.inc"
3182	#include "indicesol.inc"
3183	#include "YOMCST.inc"
3184
3185	! Parametres d'entree
3186	integer, intent(IN) :: knon, nisurf, klon
3187	real , intent(IN) :: dtime
3188	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
3189	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
3190	real, dimension(klon), intent(IN) :: ps, q1lay
3191	real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
3192	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
3193	real, dimension(klon), intent(IN) :: radsol, dif_grnd
3194	real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
3195	real, dimension(klon), intent(INOUT) :: snow, qsurf
3196
3197	! Parametres sorties
3198	real, dimension(klon), intent(OUT):: tsurf_new, evap, fluxsens, fluxlat
3199	real, dimension(klon), intent(OUT):: dflux_s, dflux_l
3200
3201	! Variables locales
3202	integer :: i
3203	real, dimension(klon) :: zx_mh, zx_nh, zx_oh
3204	real, dimension(klon) :: zx_mq, zx_nq, zx_oq
3205	real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
3206	real, dimension(klon) :: zx_sl, zx_k1
3207	real, dimension(klon) :: zx_q_0 , d_ts
3208	real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
3209	real :: bilan_f, fq_fonte
3210	REAL :: subli, fsno
3211	REAL :: qsat_new, q1_new
3212	real, parameter :: t_grnd = 271.35, t_coup = 273.15
3213	!! PB temporaire en attendant mieux pour le modele de neige
3214	REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
3215	!
3216	logical, save :: check = .false.
3217	!$OMP THREADPRIVATE(check)
3218	character (len = 20) :: modname = 'calcul_fluxs'
3219	logical, save :: fonte_neige = .false.
3220	!$OMP THREADPRIVATE(fonte_neige)
3221	real, save :: max_eau_sol = 150.0
3222	!$OMP THREADPRIVATE(max_eau_sol)
3223	character (len = 80) :: abort_message
3224	logical,save :: first = .true.,second=.false.
3225	!$OMP THREADPRIVATE(first,second)
3226
3227	if (check) write(,)'Entree ', modname,' surface = ',nisurf
3228
3229	IF (check) THEN
3230	WRITE(,)' radsol (min, max)' &
3231	& , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
3232	CALL flush(6)
3233	ENDIF
3234
3235	if (size(coastalflow) /= knon .AND. nisurf == is_ter) then
3236	write(,)'Bizarre, le nombre de points continentaux'
3237	write(,)'a change entre deux appels. J''arrete ...'
3238	abort_message='Pb run_off'
3239	call abort_gcm(modname,abort_message,1)
3240	endif
3241	!
3242	! Traitement neige et humidite du sol
3243	!
3244	!!$ WRITE(,)'test calcul_flux, surface ', nisurf
3245	!!PB test
3246	!!$ if (nisurf == is_oce) then
3247	!!$ snow = 0.
3248	!!$ qsol = max_eau_sol
3249	!!$ else
3250	!!$ where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
3251	!!$ where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
3252	!!$! snow = max(0.0, snow + (precip_snow - evap) * dtime)
3253	!!$ where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
3254	!!$ endif
3255	!!$ IF (nisurf /= is_ter) qsol = max_eau_sol
3256
3257
3258	!
3259	! Initialisation
3260	!
3261	evap = 0.
3262	fluxsens=0.
3263	fluxlat=0.
3264	dflux_s = 0.
3265	dflux_l = 0.
3266	!
3267	! zx_qs = qsat en kg/kg
3268	!
3269	DO i = 1, knon
3270	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
3271	IF (thermcep) THEN
3272	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
3273	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
3274	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
3275	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
3276	zx_qs=MIN(0.5,zx_qs)
3277	zcor=1./(1.-retv*zx_qs)
3278	zx_qs=zx_qs*zcor
3279	zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
3280	& /RLVTT / zx_pkh(i)
3281	ELSE
3282	IF (tsurf(i).LT.t_coup) THEN
3283	zx_qs = qsats(tsurf(i)) / ps(i)
3284	zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
3285	& / zx_pkh(i)
3286	ELSE
3287	zx_qs = qsatl(tsurf(i)) / ps(i)
3288	zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
3289	& / zx_pkh(i)
3290	ENDIF
3291	ENDIF
3292	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
3293	zx_qsat(i) = zx_qs
3294	zx_coef(i) = coef1lay(i) &
3295	& * (1.0+SQRT(u1lay(i)2+v1lay(i)2)) &
3296	& * p1lay(i)/(RD*t1lay(i))
3297
3298	ENDDO
3299
3300
3301	! === Calcul de la temperature de surface ===
3302	!
3303	! zx_sl = chaleur latente d'evaporation ou de sublimation
3304	!
3305	do i = 1, knon
3306	zx_sl(i) = RLVTT
3307	if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
3308	zx_k1(i) = zx_coef(i)
3309	enddo
3310
3311
3312	do i = 1, knon
3313	! Q
3314	zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
3315	zx_mq(i) = beta(i) * zx_k1(i) * &
3316	& (peqAcoef(i) - zx_qsat(i) &
3317	& + zx_dq_s_dt(i) * tsurf(i)) &
3318	& / zx_oq(i)
3319	zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
3320	& / zx_oq(i)
3321
3322	! H
3323	zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
3324	zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
3325	zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
3326
3327	! Tsurface
3328	tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
3329	& (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
3330	& + dif_grnd(i) * t_grnd * dtime)/ &
3331	& ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
3332	& zx_nh(i) + zx_sl(i) * zx_nq(i)) &
3333	& + dtime * dif_grnd(i))
3334
3335	!
3336	! Y'a-t-il fonte de neige?
3337	!
3338	! fonte_neige = (nisurf /= is_oce) .AND. &
3339	! & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
3340	! & .AND. (tsurf_new(i) >= RTT)
3341	! if (fonte_neige) tsurf_new(i) = RTT
3342	d_ts(i) = tsurf_new(i) - tsurf(i)
3343	! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
3344	! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
3345	!== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
3346	!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
3347	evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
3348	fluxlat(i) = - evap(i) * zx_sl(i)
3349	fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
3350	! Derives des flux dF/dTs (W m-2 K-1):
3351	dflux_s(i) = zx_nh(i)
3352	dflux_l(i) = (zx_sl(i) * zx_nq(i))
3353	! Nouvelle valeure de l'humidite au dessus du sol
3354	qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
3355	q1_new = peqAcoef(i) - peqBcoef(i)evap(i)dtime
3356	qsurf(i)=q1_new(1.-beta(i)) + beta(i)qsat_new
3357	!
3358	! en cas de fonte de neige
3359	!
3360	! if (fonte_neige) then
3361	! bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
3362	! & dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
3363	! & RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
3364	! bilan_f = max(0., bilan_f)
3365	! fq_fonte = bilan_f / zx_sl(i)
3366	! snow(i) = max(0., snow(i) - fq_fonte * dtime)
3367	! qsol(i) = qsol(i) + (fq_fonte * dtime)
3368	! endif
3369	!!$ if (nisurf == is_ter) &
3370	!!$ & run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
3371	!!$ qsol(i) = min(qsol(i), max_eau_sol)
3372	ENDDO
3373
3374	END SUBROUTINE calcul_fluxs
3375	!
3376	!#########################################################################
3377	!
3378	SUBROUTINE gath2cpl(champ_in, champ_out, klon, knon, iim, jmp1, knindex)
3379	use dimphy, only: liste_i,liste_j,jjphy_begin,jjphy_nb,phy_rank,phy_size
3380	implicit none
3381
3382	! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
3383	! au coupleur.
3384	!
3385	!
3386	! input:
3387	! champ_in champ sur la grille gathere
3388	! knon nombre de points dans le domaine a traiter
3389	! knindex index des points de la surface a traiter
3390	! klon taille de la grille
3391	! iim,jjm dimension de la grille 2D
3392	!
3393	! output:
3394	! champ_out champ sur la grille 2D
3395	!
3396	! input
3397	integer :: klon, knon, iim, jmp1
3398	real, dimension(klon) :: champ_in
3399	integer, dimension(klon) :: knindex
3400	! output
3401	real, dimension(iim,jmp1) :: champ_out
3402	! local
3403	integer :: i, ig, j
3404	real, dimension(klon) :: tamp
3405
3406	tamp = 0.
3407	do i = 1, knon
3408	ig = knindex(i)
3409	tamp(ig) = champ_in(i)
3410	enddo
3411
3412	!ym ig = 1
3413	!ym champ_out(:,1) = tamp(ig)
3414	!ym do j = 2, jjm
3415	!ym do i = 1, iim
3416	!ym ig = ig + 1
3417	!ym champ_out(i,j) = tamp(ig)
3418	!ym enddo
3419	!ym enddo
3420	!ym ig = ig + 1
3421	!ym champ_out(:,jjm+1) = tamp(ig)
3422
3423	do ig=1,klon
3424	i=liste_i(ig)
3425	j=liste_j(ig)-jjphy_begin+1
3426	champ_out(i,j)=tamp(ig)
3427	enddo
3428
3429	if (phy_rank==0) champ_out(:,1)=tamp(1)
3430	if (phy_rank==phy_size-1) champ_out(:,jjphy_nb)=tamp(klon)
3431
3432	END SUBROUTINE gath2cpl
3433	!
3434	!#########################################################################
3435	!
3436	SUBROUTINE cpl2gath(champ_in, champ_out, klon, knon, iim, jmp1, knindex)
3437	use dimphy, only : liste_i, liste_j, jjphy_begin
3438	implicit none
3439	! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
3440	! au coupleur.
3441	!
3442	!
3443	! input:
3444	! champ_in champ sur la grille gathere
3445	! knon nombre de points dans le domaine a traiter
3446	! knindex index des points de la surface a traiter
3447	! klon taille de la grille
3448	! iim,jjm dimension de la grille 2D
3449	!
3450	! output:
3451	! champ_out champ sur la grille 2D
3452	!
3453	! input
3454	integer :: klon, knon, iim, jmp1
3455	real, dimension(iim,jmp1) :: champ_in
3456	integer, dimension(klon) :: knindex
3457	! output
3458	real, dimension(klon) :: champ_out
3459	! local
3460	integer :: i, ig, j
3461	real, dimension(klon) :: tamp
3462	logical ,save :: check = .false.
3463
3464	!ym ig = 1
3465	!ym tamp(ig) = champ_in(1,1)
3466	!ym do j = 2, jjm
3467	!ym do i = 1, iim
3468	!ym ig = ig + 1
3469	!ym tamp(ig) = champ_in(i,j)
3470	!ym enddo
3471	!ym enddo
3472	!ym ig = ig + 1
3473	!ym tamp(ig) = champ_in(1,jjm+1)
3474
3475	do ig=1,klon
3476	i=liste_i(ig)
3477	j=liste_j(ig)-jjphy_begin+1
3478	tamp(ig)=champ_in(i,j)
3479	enddo
3480
3481	do i = 1, knon
3482	ig = knindex(i)
3483	champ_out(i) = tamp(ig)
3484	enddo
3485
3486	END SUBROUTINE cpl2gath
3487	!
3488	!#########################################################################
3489	!
3490	SUBROUTINE albsno(klon, knon,dtime,agesno,alb_neig_grid, precip_snow)
3491	IMPLICIT none
3492
3493	INTEGER :: klon, knon
3494	INTEGER, PARAMETER :: nvm = 8
3495	REAL :: dtime
3496	REAL, dimension(klon,nvm) :: veget
3497	REAL, DIMENSION(klon) :: alb_neig_grid, agesno, precip_snow
3498
3499	INTEGER :: i, nv
3500
3501	REAL, DIMENSION(nvm),SAVE :: init, decay
3502	!$OMP THREADPRIVATE(init, decay)
3503	REAL :: as
3504	DATA init /0.55, 0.14, 0.18, 0.29, 0.15, 0.15, 0.14, 0./
3505	DATA decay/0.30, 0.67, 0.63, 0.45, 0.40, 0.14, 0.06, 1./
3506
3507	veget = 0.
3508	veget(:,1) = 1. ! desert partout
3509	DO i = 1, knon
3510	alb_neig_grid(i) = 0.0
3511	ENDDO
3512	DO nv = 1, nvm
3513	DO i = 1, knon
3514	as = init(nv)+decay(nv)*EXP(-agesno(i)/5.)
3515	alb_neig_grid(i) = alb_neig_grid(i) + veget(i,nv)*as
3516	ENDDO
3517	ENDDO
3518	!
3519	!! modilation en fonction de l'age de la neige
3520	!
3521	DO i = 1, knon
3522	agesno(i) = (agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
3523	& * EXP(-1.MAX(0.0,precip_snow(i))dtime/0.3)
3524	agesno(i) = MAX(agesno(i),0.0)
3525	ENDDO
3526
3527	END SUBROUTINE albsno
3528	!
3529	!#########################################################################
3530	!
3531
3532	SUBROUTINE fonte_neige( klon, knon, nisurf, dtime, &
3533	& tsurf, p1lay, cal, beta, coef1lay, ps, &
3534	& precip_rain, precip_snow, snow, qsol, &
3535	& radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
3536	& petAcoef, peqAcoef, petBcoef, peqBcoef, &
3537	& tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
3538	& fqcalving,fqfonte,ffonte,run_off_lic_0)
3539
3540	! Routine de traitement de la fonte de la neige dans le cas du traitement
3541	! de sol simplifi�!
3542	! LF 03/2001
3543	! input:
3544	! knon nombre de points a traiter
3545	! nisurf surface a traiter
3546	! tsurf temperature de surface
3547	! p1lay pression 1er niveau (milieu de couche)
3548	! cal capacite calorifique du sol
3549	! beta evap reelle
3550	! coef1lay coefficient d'echange
3551	! ps pression au sol
3552	! precip_rain precipitations liquides
3553	! precip_snow precipitations solides
3554	! snow champs hauteur de neige
3555	! qsol hauteur d'eau contenu dans le sol
3556	! runoff runoff en cas de trop plein
3557	! petAcoef coeff. A de la resolution de la CL pour t
3558	! peqAcoef coeff. A de la resolution de la CL pour q
3559	! petBcoef coeff. B de la resolution de la CL pour t
3560	! peqBcoef coeff. B de la resolution de la CL pour q
3561	! radsol rayonnement net aus sol (LW + SW)
3562	! dif_grnd coeff. diffusion vers le sol profond
3563	!
3564	! output:
3565	! tsurf_new temperature au sol
3566	! fluxsens flux de chaleur sensible
3567	! fluxlat flux de chaleur latente
3568	! dflux_s derivee du flux de chaleur sensible / Ts
3569	! dflux_l derivee du flux de chaleur latente / Ts
3570	! in/out:
3571	! run_off_lic_0 run off glacier du pas de temps pr�edent
3572	!
3573
3574	#include "YOETHF.inc"
3575	!rv#include "FCTTRE.inc"
3576	#include "indicesol.inc"
3577	!IM cf JLD
3578	#include "YOMCST.inc"
3579
3580	! Parametres d'entree
3581	integer, intent(IN) :: knon, nisurf, klon
3582	real , intent(IN) :: dtime
3583	real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
3584	real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
3585	real, dimension(klon), intent(IN) :: ps, q1lay
3586	real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
3587	real, dimension(klon), intent(IN) :: precip_rain, precip_snow
3588	real, dimension(klon), intent(IN) :: radsol, dif_grnd
3589	real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
3590	real, dimension(klon), intent(INOUT) :: snow, qsol
3591
3592	! Parametres sorties
3593	real, dimension(klon), intent(INOUT):: tsurf_new, evap, fluxsens, fluxlat
3594	real, dimension(klon), intent(INOUT):: dflux_s, dflux_l
3595	! Flux thermique utiliser pour fondre la neige
3596	real, dimension(klon), intent(INOUT):: ffonte
3597	! Flux d'eau "perdu" par la surface et necessaire pour que limiter la
3598	! hauteur de neige, en kg/m2/s. Et flux d'eau de fonte de la calotte.
3599	REAL, DIMENSION(klon), INTENT(INOUT):: fqcalving, fqfonte
3600	real, dimension(klon), intent(INOUT):: run_off_lic_0
3601	! Variables locales
3602	! Masse maximum de neige (kg/m2). Au dessus de ce seuil, la neige
3603	! en exces "s'ecoule" (calving)
3604	! real, parameter :: snow_max=1.
3605	!IM cf JLD/GK
3606	real, parameter :: snow_max=3000.
3607	integer :: i
3608	real, dimension(klon) :: zx_mh, zx_nh, zx_oh
3609	real, dimension(klon) :: zx_mq, zx_nq, zx_oq
3610	real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
3611	real, dimension(klon) :: zx_sl, zx_k1
3612	real, dimension(klon) :: zx_q_0 , d_ts
3613	real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
3614	real :: bilan_f, fq_fonte
3615	REAL :: subli, fsno
3616	REAL, DIMENSION(klon) :: bil_eau_s, snow_evap
3617	real, parameter :: t_grnd = 271.35, t_coup = 273.15
3618	!! PB temporaire en attendant mieux pour le modele de neige
3619	! REAL, parameter :: chasno = RLMLT/(2.3867E+06*0.15)
3620	REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
3621	!IM cf JLD/ GKtest
3622	REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)
3623	! fin GKtest
3624	!
3625	logical, save :: check = .FALSE.
3626	!$OMP THREADPRIVATE(check)
3627	character (len = 20) :: modname = 'fonte_neige'
3628	logical, save :: neige_fond = .false.
3629	!$OMP THREADPRIVATE(neige_fond)
3630	real, save :: max_eau_sol = 150.0
3631	!$OMP THREADPRIVATE(max_eau_sol)
3632	character (len = 80) :: abort_message
3633	logical,save :: first = .true.,second=.false.
3634	!$OMP THREADPRIVATE(first,second)
3635	real :: coeff_rel
3636	#include "FCTTRE.inc"
3637
3638
3639	if (check) write(,)'Entree ', modname,' surface = ',nisurf
3640
3641	! Initialisations
3642	coeff_rel = dtime/(tau_calv * rday)
3643	bil_eau_s(:) = 0.
3644	DO i = 1, knon
3645	zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
3646	IF (thermcep) THEN
3647	zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
3648	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
3649	zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
3650	zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
3651	zx_qs=MIN(0.5,zx_qs)
3652	zcor=1./(1.-retv*zx_qs)
3653	zx_qs=zx_qs*zcor
3654	zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
3655	& /RLVTT / zx_pkh(i)
3656	ELSE
3657	IF (tsurf(i).LT.t_coup) THEN
3658	zx_qs = qsats(tsurf(i)) / ps(i)
3659	zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
3660	& / zx_pkh(i)
3661	ELSE
3662	zx_qs = qsatl(tsurf(i)) / ps(i)
3663	zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
3664	& / zx_pkh(i)
3665	ENDIF
3666	ENDIF
3667	zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
3668	zx_qsat(i) = zx_qs
3669	zx_coef(i) = coef1lay(i) &
3670	& * (1.0+SQRT(u1lay(i)2+v1lay(i)2)) &
3671	& * p1lay(i)/(RD*t1lay(i))
3672	ENDDO
3673
3674
3675	! === Calcul de la temperature de surface ===
3676	!
3677	! zx_sl = chaleur latente d'evaporation ou de sublimation
3678	!
3679	do i = 1, knon
3680	zx_sl(i) = RLVTT
3681	if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
3682	zx_k1(i) = zx_coef(i)
3683	enddo
3684
3685
3686	do i = 1, knon
3687	! Q
3688	zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
3689	zx_mq(i) = beta(i) * zx_k1(i) * &
3690	& (peqAcoef(i) - zx_qsat(i) &
3691	& + zx_dq_s_dt(i) * tsurf(i)) &
3692	& / zx_oq(i)
3693	zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
3694	& / zx_oq(i)
3695
3696	! H
3697	zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
3698	zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
3699	zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
3700	enddo
3701
3702
3703	WHERE (precip_snow > 0.) snow = snow + (precip_snow * dtime)
3704	snow_evap = 0.
3705	WHERE (evap > 0. )
3706	snow_evap = MIN (snow / dtime, evap)
3707	snow = snow - snow_evap * dtime
3708	snow = MAX(0.0, snow)
3709	end where
3710
3711	! bil_eau_s = bil_eau_s + (precip_rain * dtime) - (evap - snow_evap) * dtime
3712	bil_eau_s = (precip_rain * dtime) - (evap - snow_evap) * dtime
3713
3714	!
3715	! Y'a-t-il fonte de neige?
3716	!
3717	ffonte=0.
3718	do i = 1, knon
3719	neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
3720	& .AND. tsurf_new(i) >= RTT)
3721	if (neige_fond) then
3722	fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i))
3723	ffonte(i) = fq_fonte * RLMLT/dtime
3724	fqfonte(i) = fq_fonte/dtime
3725	snow(i) = max(0., snow(i) - fq_fonte)
3726	bil_eau_s(i) = bil_eau_s(i) + fq_fonte
3727	tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno
3728	!IM cf JLD OK
3729	!IM cf JLD/ GKtest fonte aussi pour la glace
3730	IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
3731	fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,0.0)
3732	ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime
3733	IF ( ok_lic_melt ) THEN
3734	fqfonte(i) = fqfonte(i) + fq_fonte/dtime
3735	bil_eau_s(i) = bil_eau_s(i) + fq_fonte
3736	ENDIF
3737	tsurf_new(i) = RTT
3738	ENDIF
3739	d_ts(i) = tsurf_new(i) - tsurf(i)
3740	endif
3741	!
3742	! s'il y a une hauteur trop importante de neige, elle s'coule
3743	fqcalving(i) = max(0., snow(i) - snow_max)/dtime
3744	snow(i)=min(snow(i),snow_max)
3745	!
3746	IF (nisurf == is_ter) then
3747	qsol(i) = qsol(i) + bil_eau_s(i)
3748	run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)
3749	qsol(i) = MIN(qsol(i), max_eau_sol)
3750	else if (nisurf == is_lic) then
3751	run_off_lic(i) = (coeff_rel * fqcalving(i)) + &
3752	& (1. - coeff_rel) * run_off_lic_0(i)
3753	run_off_lic_0(i) = run_off_lic(i)
3754	run_off_lic(i) = run_off_lic(i) + fqfonte(i)/dtime
3755	endif
3756	enddo
3757
3758	END SUBROUTINE fonte_neige
3759	!
3760	!#########################################################################
3761	!
3762	END MODULE interface_surf

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

Download in other formats:

Original Format