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

Last change on this file was 564, checked in by lmdzadmin, 20 years ago

Encore quelques initialisations pour le portage AC
LF

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