Context Navigation

interface_surf.F90 @ 279

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

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

Download in other formats:

Original Format