Context Navigation

interface_surf.F90 @ 181

Last change on this file since 181 was 181, checked in by lmdzadmin, 24 years ago
Traitement de la fonte de la neige dans le cas du modele de sol simplifie est sorti de calcul_flux LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 82.9 KB

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

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

Download in other formats:

Original Format