Context Navigation

interface_surf.F90 @ 1566

Last change on this file since 1566 was 634, checked in by Laurent Fairhead, 20 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 99.4 KB

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

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

Download in other formats:

Original Format