source: LMDZ.3.3/branches/rel-LF/libf/phylmd/interface_surf.F90 @ 436

Last change on this file since 436 was 398, checked in by lmdzadmin, 22 years ago

Réglage decalage runoff entre orchidee et coupleur OM
LF

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