Context Navigation

interface_surf.F90 @ 281

Last change on this file since 281 was 281, checked in by lmdzadmin, 23 years ago
Passage des deux albedos de surface (vis et nir) Modif dans interfsol sur les indicages des variables passees a ORCHIDEE pour ne plus avoir de decalage dans les sorties ORCHIDEE LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 88.3 KB

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

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

Download in other formats:

Original Format