source: LMDZ4/branches/V3_test/libf/phylmd/interface_surf.F90

Last change on this file was 752, checked in by lsce, 18 years ago

ACo + AC : efface un print inutile et encombrant

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 123.0 KB
Line 
1!
2! $Header$
3!
4
5  MODULE interface_surf
6
7! Ce module regroupe toutes les routines gerant l'interface entre le modele
8! atmospherique et les modeles de surface (sols continentaux, oceans, glaces)
9! Les routines sont les suivantes:
10!
11!   interfsurf_*: routines d'aiguillage vers les interfaces avec les
12!                 differents modeles de surface
13!   interfsol\
14!             > routines d'interface proprement dite
15!   interfoce/
16!
17!   interfstart: routine d'initialisation et de lecture de l'etat initial
18!                "interface"
19!   interffin  : routine d'ecriture de l'etat de redemmarage de l'interface
20!
21!
22! L. Fairhead, LMD, 02/2000
23
24!ym  USE ioipsl
25
26  IMPLICIT none
27
28  PRIVATE
29  PUBLIC :: interfsurf,interfsurf_hq, gath2cpl
30
31  INTERFACE interfsurf
32    module procedure interfsurf_hq, interfsurf_vent
33  END INTERFACE
34
35  INTERFACE interfoce
36    module procedure interfoce_cpl, interfoce_slab, interfoce_lim
37  END INTERFACE
38
39#include "YOMCST.inc"
40#include "indicesol.inc"
41!IM
42#include "clesphys.inc"
43
44! run_off      ruissellement total
45  REAL, ALLOCATABLE, DIMENSION(:),SAVE    :: run_off, run_off_lic
46!$OMP THREADPRIVATE(run_off, run_off_lic)
47  real, allocatable, dimension(:),save    :: coastalflow, riverflow
48!$OMP THREADPRIVATE(coastalflow, riverflow)
49!!$PB
50  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: tmp_rriv, tmp_rcoa,tmp_rlic
51!$OMP THREADPRIVATE(tmp_rriv, tmp_rcoa,tmp_rlic)
52!! pour simuler la fonte des glaciers antarctiques
53  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: coeff_iceberg
54!$OMP THREADPRIVATE(coeff_iceberg)
55  real, save                              :: surf_maille
56!$OMP THREADPRIVATE(surf_maille)
57  real, save                              :: cte_flux_iceberg = 6.3e7
58!$OMP THREADPRIVATE(cte_flux_iceberg) 
59  integer, save                           :: num_antarctic = 1
60!$OMP THREADPRIVATE(num_antarctic)
61  REAL, save                              :: tau_calv
62!$OMP THREADPRIVATE(tau_calv)
63!!$
64  CONTAINS
65!
66!############################################################################
67!
68  SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, &
69      & klon, iim, jjm, nisurf, knon, knindex, pctsrf, &
70      & rlon, rlat, cufi, cvfi,&
71      & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, qsol,&
72      & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
73      & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
74      & precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &
75      & fder, taux, tauy, &
76! -- LOOP
77      & windsp, &
78! -- LOOP
79      & rugos, rugoro, &
80      & albedo, snow, qsurf, &
81      & tsurf, p1lay, ps, radsol, &
82      & ocean, npas, nexca, zmasq, &
83      & evap, fluxsens, fluxlat, dflux_l, dflux_s, &             
84      & tsol_rad, tsurf_new, alb_new, alblw, emis_new, &
85      & z0_new, pctsrf_new, agesno,fqcalving,fqfonte,ffonte, run_off_lic_0,&
86!IM "slab" ocean
87      & flux_o, flux_g, tslab, seaice)
88
89
90   USE dimphy,only : monocpu,jjphy_nb,omp_rank
91! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general
92! (sols continentaux, oceans, glaces) pour les fluxs de chaleur et d'humidite.
93! En pratique l'interface se fait entre la couche limite du modele
94! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
95!
96!
97! L.Fairhead 02/2000
98!
99! input:
100!   itime        numero du pas de temps
101!   klon         nombre total de points de grille
102!   iim, jjm     nbres de pts de grille
103!   dtime        pas de temps de la physique (en s)
104!   date0        jour initial
105!   jour         jour dans l'annee en cours,
106!   rmu0         cosinus de l'angle solaire zenithal
107!   nexca        pas de temps couplage
108!   nisurf       index de la surface a traiter (1 = sol continental)
109!   knon         nombre de points de la surface a traiter
110!   knindex      index des points de la surface a traiter
111!   pctsrf       tableau des pourcentages de surface de chaque maille
112!   rlon         longitudes
113!   rlat         latitudes
114!   cufi,cvfi    resolution des mailles en x et y (m)
115!   debut        logical: 1er appel a la physique
116!   lafin        logical: dernier appel a la physique
117!   ok_veget     logical: appel ou non au schema de surface continental
118!                     (si false calcul simplifie des fluxs sur les continents)
119!   zlev         hauteur de la premiere couche
120!   u1_lay       vitesse u 1ere couche
121!   v1_lay       vitesse v 1ere couche
122!   temp_air     temperature de l'air 1ere couche
123!   spechum      humidite specifique 1ere couche
124!   epot_air     temp potentielle de l'air
125!   ccanopy      concentration CO2 canopee
126!   tq_cdrag     cdrag
127!   petAcoef     coeff. A de la resolution de la CL pour t
128!   peqAcoef     coeff. A de la resolution de la CL pour q
129!   petBcoef     coeff. B de la resolution de la CL pour t
130!   peqBcoef     coeff. B de la resolution de la CL pour q
131!   precip_rain  precipitation liquide
132!   precip_snow  precipitation solide
133!   sollw        flux IR net a la surface
134!   sollwdown    flux IR descendant a la surface
135!   swnet        flux solaire net
136!   swdown       flux solaire entrant a la surface
137!   albedo       albedo de la surface
138!   tsurf        temperature de surface
139!   tslab        temperature slab ocean
140!   pctsrf_slab  pourcentages (0-1) des sous-surfaces dans le slab
141!   tmp_pctsrf_slab = pctsrf_slab
142!   p1lay        pression 1er niveau (milieu de couche)
143!   ps           pression au sol
144!   radsol       rayonnement net aus sol (LW + SW)
145!   ocean        type d'ocean utilise (force, slab, couple)
146!   fder         derivee des flux (pour le couplage)
147!   taux, tauy   tension de vents
148! -- LOOP
149!   windsp       module du vent a 10m
150! -- LOOP
151!   rugos        rugosite
152!   zmasq        masque terre/ocean
153!   rugoro       rugosite orographique
154!   run_off_lic_0 runoff glacier du pas de temps precedent
155!
156! output:
157!   evap         evaporation totale
158!   fluxsens     flux de chaleur sensible
159!   fluxlat      flux de chaleur latente
160!   tsol_rad     
161!   tsurf_new    temperature au sol
162!   alb_new      albedo
163!   emis_new     emissivite
164!   z0_new       surface roughness
165!   pctsrf_new   nouvelle repartition des surfaces
166
167#include "iniprint.h"
168
169
170! Parametres d'entree
171  integer, intent(IN) :: itime
172  integer, intent(IN) :: iim, jjm
173  integer, intent(IN) :: klon
174  real, intent(IN) :: dtime
175  real, intent(IN) :: date0
176  integer, intent(IN) :: jour
177  real, intent(IN)    :: rmu0(klon)
178  integer, intent(IN) :: nisurf
179  integer, intent(IN) :: knon
180  integer, dimension(klon), intent(in) :: knindex
181  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
182  logical, intent(IN) :: debut, lafin, ok_veget
183  real, dimension(klon), intent(IN) :: rlon, rlat
184  real, dimension(klon), intent(IN) :: cufi, cvfi
185  real, dimension(klon), intent(INOUT) :: tq_cdrag
186  real, dimension(klon), intent(IN) :: zlev
187  real, dimension(klon), intent(IN) :: u1_lay, v1_lay
188  real, dimension(klon), intent(IN) :: temp_air, spechum
189  real, dimension(klon), intent(IN) :: epot_air, ccanopy
190  real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
191  real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
192  real, dimension(klon), intent(IN) :: precip_rain, precip_snow
193  real, dimension(klon), intent(IN) :: sollw, sollwdown, swnet, swdown
194  real, dimension(klon), intent(IN) :: ps, albedo
195  real, dimension(klon), intent(IN) :: tsurf, p1lay
196!IM: "slab" ocean
197  real :: amn, amx
198  real, dimension(klon), intent(INOUT) :: tslab
199  real, allocatable, dimension(:), save :: tmp_tslab
200!$OMP THREADPRIVATE(tmp_tslab)
201  real, dimension(klon), intent(OUT) :: flux_o, flux_g
202  real, dimension(klon), intent(INOUT)        :: seaice ! glace de mer (kg/m2)
203  real, dimension(klon)                       :: siceh  ! hauteur glace de mer (m)
204  REAL, DIMENSION(klon), INTENT(INOUT) :: radsol,fder
205 
206!  real, dimension(klon), intent(IN) :: zmasq
207  real, dimension(klon), intent(IN) :: zmasq
208  real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro
209! -- LOOP
210  real, dimension(klon), intent(IN) :: windsp
211! -- LOOP
212  character (len = 6)  :: ocean
213  integer              :: npas, nexca ! nombre et pas de temps couplage
214  real, dimension(klon), intent(INOUT) :: evap, snow, qsurf
215!! PB ajout pour soil
216  logical          :: soil_model
217  integer          :: nsoilmx
218  REAL, DIMENSION(klon, nsoilmx) :: tsoil
219  REAL, dimension(klon), intent(INOUT) :: qsol
220  REAL, dimension(klon)          :: soilcap
221  REAL, dimension(klon)          :: soilflux
222! Parametres de sortie
223  real, dimension(klon), intent(OUT):: fluxsens, fluxlat
224  real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new
225  real, dimension(klon), intent(OUT):: alblw
226  real, dimension(klon), intent(OUT):: emis_new, z0_new
227  real, dimension(klon), intent(OUT):: dflux_l, dflux_s
228  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
229  real, dimension(klon), intent(INOUT):: agesno
230  real, dimension(klon), intent(INOUT):: run_off_lic_0
231
232! Flux thermique utiliser pour fondre la neige
233!jld a rajouter   real, dimension(klon), intent(INOUT):: ffonte
234  real, dimension(klon), intent(INOUT):: ffonte
235! Flux d'eau "perdue" par la surface et necessaire pour que limiter la
236! hauteur de neige, en kg/m2/s. Et quantite d'eau de fonte de la calotte.
237!jld a rajouter   real, dimension(klon), intent(INOUT):: fqcalving, fqfonte
238  REAL, DIMENSION(klon), INTENT(INOUT):: fqcalving, fqfonte
239!IM: "slab" ocean
240  real, dimension(klon) :: new_dif_grnd
241!IM: "slab" ocean - Local
242  real, parameter :: t_grnd=271.35
243  real, dimension(klon) :: zx_sl
244  integer i
245  real, allocatable, dimension(:), save :: tmp_flux_o, tmp_flux_g
246!$OMP THREADPRIVATE(tmp_flux_o, tmp_flux_g) 
247  real, allocatable, dimension(:), save :: tmp_radsol
248!$OMP THREADPRIVATE(tmp_radsol)
249  real, allocatable, dimension(:,:), save :: tmp_pctsrf_slab
250!$OMP THREADPRIVATE(tmp_pctsrf_slab)
251  real, allocatable, dimension(:), save :: tmp_seaice
252!$OMP THREADPRIVATE(tmp_seaice)
253! Local
254  character (len = 20),save :: modname = 'interfsurf_hq'
255!$OMP THREADPRIVATE(modname)
256  character (len = 80) :: abort_message
257  logical, save        :: first_call = .true.
258!$OMP THREADPRIVATE(first_call) 
259  integer, save        :: error
260!$OMP THREADPRIVATE(error) 
261  integer              :: ii, index
262  logical,save              :: check = .true.
263!$OMP THREADPRIVATE(check) 
264  real, dimension(klon):: cal, beta, dif_grnd, capsol
265!!$PB  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.
266  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.
267  real, parameter      :: calsno=1./(2.3867e+06*.15)
268  real, dimension(klon):: alb_ice
269  real, dimension(klon):: tsurf_temp
270  real, dimension(klon):: qsurf_new
271!!  real, allocatable, dimension(:), save :: alb_neig_grid
272  real, dimension(klon):: alb_neig, alb_eau
273  real, DIMENSION(klon):: zfra
274  logical              :: cumul = .false.
275  INTEGER,dimension(1) :: iloc
276  INTEGER                 :: isize
277  real, dimension(klon):: fder_prev
278  REAL, dimension(klon) :: bidule
279  real, dimension(klon) :: ps_tmp,p1lay_tmp
280  INTEGER :: j
281!
282!IM ?? quelques variables pour netcdf
283#include "netcdf.inc"
284
285  if (check) write(*,*) 'Entree ', modname
286!
287! On doit commencer par appeler les schemas de surfaces continentales
288! car l'ocean a besoin du ruissellement qui est y calcule
289!
290  if (first_call) then
291    call conf_interface(tau_calv)
292    if (nisurf /= is_ter .and. klon > 1) then
293      write(*,*)' *** Warning ***'
294      write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter
295      write(*,*)'or on doit commencer par les surfaces continentales'
296      abort_message='voir ci-dessus'
297      call abort_gcm(modname,abort_message,1)
298    endif
299    if (ocean /= 'slab  ' .and. ocean /= 'force ' .and. ocean /= 'couple') then
300      write(*,*)' *** Warning ***'
301      write(*,*)'Option couplage pour l''ocean = ', ocean
302      abort_message='option pour l''ocean non valable'
303      call abort_gcm(modname,abort_message,1)
304    endif
305    if ( is_oce > is_sic ) then
306      write(*,*)' *** Warning ***'
307      write(*,*)' Pour des raisons de sequencement dans le code'
308      write(*,*)' l''ocean doit etre traite avant la banquise'
309      write(*,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic
310      abort_message='voir ci-dessus'
311      call abort_gcm(modname,abort_message,1)
312    endif
313!    allocate(alb_neig_grid(klon), stat = error)
314!    if (error /= 0) then
315!      abort_message='Pb allocation alb_neig_grid'
316!      call abort_gcm(modname,abort_message,1)
317!    endif
318  endif
319  first_call = .false.
320 
321! Initialisations diverses
322!
323!!$  cal=0.; beta=1.; dif_grnd=0.; capsol=0.
324!!$  alb_new = 0.; z0_new = 0.; alb_neig = 0.0
325!!$! PB
326!!$  tsurf_new = 0.
327
328!IM cf JLD
329  ffonte(1:knon)=0.
330  fqcalving(1:knon)=0.
331  fqfonte  (1:knon)=0.
332
333
334  cal = 999999. ; beta = 999999. ; dif_grnd = 999999. ; capsol = 999999.
335  alb_new = 999999. ; z0_new = 999999. ; alb_neig = 999999.
336  tsurf_new = 999999.
337  alblw = 999999.
338
339!IM: "slab" ocean; initialisations
340  flux_o = 0.
341  flux_g = 0.
342!
343    if (.not. allocated(tmp_flux_o)) then
344      allocate(tmp_flux_o(klon), stat = error)
345      DO i=1, knon
346       tmp_flux_o(knindex(i))=flux_o(i)
347      ENDDO
348      if (error /= 0) then
349        abort_message='Pb allocation tmp_flux_o'
350        call abort_gcm(modname,abort_message,1)
351      endif
352    endif
353    if (.not. allocated(tmp_flux_g)) then 
354      allocate(tmp_flux_g(klon), stat = error)
355      DO i=1, knon
356       tmp_flux_g(knindex(i))=flux_g(i)
357      ENDDO
358      if (error /= 0) then
359        abort_message='Pb allocation tmp_flux_g'
360        call abort_gcm(modname,abort_message,1)
361      endif
362    endif
363    if (.not. allocated(tmp_radsol)) then 
364      allocate(tmp_radsol(klon), stat = error)
365      if (error /= 0) then
366        abort_message='Pb allocation tmp_radsol'
367        call abort_gcm(modname,abort_message,1)
368      endif
369    endif
370    if (.not. allocated(tmp_pctsrf_slab)) then
371      allocate(tmp_pctsrf_slab(klon,nbsrf), stat = error)
372      if (error /= 0) then
373        abort_message='Pb allocation tmp_pctsrf_slab'
374        call abort_gcm(modname,abort_message,1)
375      endif
376    DO i=1, klon
377     tmp_pctsrf_slab(i,1:nbsrf)=pctsrf(i,1:nbsrf)
378    ENDDO
379    endif
380!
381    if (.not. allocated(tmp_seaice)) then
382      allocate(tmp_seaice(klon), stat = error)
383      if (error /= 0) then
384        abort_message='Pb allocation tmp_seaice'
385        call abort_gcm(modname,abort_message,1)
386      endif
387    IF(check) THEN
388     PRINT*,'allocation tmp_seaice nisurf itime',nisurf, itime
389    ENDIF
390    endif
391!
392    if (.not. allocated(tmp_tslab)) then
393      allocate(tmp_tslab(klon), stat = error)
394      if (error /= 0) then
395        abort_message='Pb allocation tmp_tslab'
396        call abort_gcm(modname,abort_message,1)
397      endif
398    endif
399    DO i=1, klon
400     tmp_tslab(i)=tslab(i)
401    ENDDO
402!
403! Aiguillage vers les differents schemas de surface
404
405  if (nisurf == is_ter) then
406!
407! Surface "terre" appel a l'interface avec les sols continentaux
408!
409! allocation du run-off
410    if (.not. allocated(coastalflow)) then
411      allocate(coastalflow(knon), stat = error)
412      if (error /= 0) then
413        abort_message='Pb allocation coastalflow'
414        call abort_gcm(modname,abort_message,1)
415      endif
416      allocate(riverflow(knon), stat = error)
417      if (error /= 0) then
418        abort_message='Pb allocation riverflow'
419        call abort_gcm(modname,abort_message,1)
420      endif
421      allocate(run_off(knon), stat = error)
422      if (error /= 0) then
423        abort_message='Pb allocation run_off'
424        call abort_gcm(modname,abort_message,1)
425      endif
426!cym     
427      run_off=0.0
428!cym
429
430!!$PB
431      ALLOCATE (tmp_rriv(iim,jjphy_nb), stat=error)
432      if (error /= 0) then
433        abort_message='Pb allocation tmp_rriv'
434        call abort_gcm(modname,abort_message,1)
435      endif
436      tmp_rriv=0.
437      ALLOCATE (tmp_rcoa(iim,jjphy_nb), stat=error)
438      if (error /= 0) then
439        abort_message='Pb allocation tmp_rcoa'
440        call abort_gcm(modname,abort_message,1)
441      endif
442      tmp_rcoa=0.
443!ym      ALLOCATE (tmp_rlic(iim,jjm+1), stat=error)
444      ALLOCATE (tmp_rlic(iim,jjphy_nb), stat=error)
445      if (error /= 0) then
446        abort_message='Pb allocation tmp_rlic'
447        call abort_gcm(modname,abort_message,1)
448      endif
449      tmp_rlic=0.
450!!$
451    else if (size(coastalflow) /= knon) then
452      write(*,*)'Bizarre, le nombre de points continentaux'
453      write(*,*)'a change entre deux appels. J''arrete ...'
454      abort_message='voir ci-dessus'
455      call abort_gcm(modname,abort_message,1)
456    endif
457    coastalflow = 0.
458    riverflow = 0.
459!
460! Calcul age de la neige
461!
462!!$ PB ATTENTION changement ordre des appels
463!!$    CALL albsno(klon,agesno,alb_neig_grid) 
464
465    if (.not. ok_veget) then
466!
467! calcul albedo: lecture albedo fichier CL puis ajout albedo neige
468!
469       call interfsur_lim(itime, dtime, jour, &
470     & klon, nisurf, knon, knindex, debut,  &
471     & alb_new, z0_new)
472
473! calcul snow et qsurf, hydrol adapte !
474       CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
475
476       IF (soil_model) THEN
477           CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)
478           cal(1:knon) = RCPD / soilcap(1:knon)
479           radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
480       ELSE
481           cal = RCPD * capsol
482!!$      cal = capsol
483       ENDIF
484       CALL calcul_fluxs( klon, knon, nisurf, dtime, &
485     &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
486     &   precip_rain, precip_snow, snow, qsurf,  &
487     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
488     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
489     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
490
491       CALL fonte_neige( klon, knon, nisurf, dtime, &
492     &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
493     &   precip_rain, precip_snow, snow, qsol,  &
494     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
495     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
496     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
497     &   fqcalving,fqfonte,ffonte, run_off_lic_0)
498
499
500     call albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
501     where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
502     zfra(1:knon) = max(0.0,min(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
503     alb_new(1 : knon)  = alb_neig(1 : knon) *zfra(1:knon) + &
504    &                     alb_new(1 : knon)*(1.0-zfra(1:knon))
505     z0_new = sqrt(z0_new**2+rugoro**2)
506     alblw(1 : knon) = alb_new(1 : knon)
507
508    else
509!!      CALL albsno(klon,agesno,alb_neig_grid) 
510!
511!  appel a sechiba
512!
513#ifdef CPP_VEGET
514      p1lay_tmp(1:knon)=p1lay(1:knon)/100.
515      ps_tmp(1:knon)=ps(1:knon)/100.
516 
517       call interfsol(itime, klon, dtime, date0, nisurf, knon, &
518     &  knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
519     &  debut, lafin, ok_veget, &
520     &  zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
521     &  tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
522     &  precip_rain, precip_snow, sollwdown, swnet, swdown, &
523     &  tsurf, p1lay_tmp, ps_tmp, radsol, &
524     &  evap, fluxsens, fluxlat, &             
525     &  tsol_rad, tsurf_new, alb_new, alblw, &
526     &  emis_new, z0_new, dflux_l, dflux_s, qsurf_new)
527
528
529! ajout de la contribution du relief
530
531      z0_new = SQRT(z0_new**2+rugoro**2)
532!
533! mise a jour de l'humidite saturante calculee par ORCHIDEE
534      qsurf(1:knon) = qsurf_new(1:knon)
535#else
536      abort_message='Pb de coherence: ok_veget = .t. mais CPP_VEGET = .f.'
537      call abort_gcm(modname,abort_message,1)
538#endif
539
540    endif   
541!
542! Remplissage des pourcentages de surface
543!
544    pctsrf_new(:,nisurf) = pctsrf(:,nisurf)
545
546  else if (nisurf == is_oce) then
547
548    if (check) write(*,*)'ocean, nisurf = ',nisurf,'knon=',knon
549!
550! Surface "ocean" appel a l'interface avec l'ocean
551!
552    if (ocean == 'couple') then
553      if (nexca == 0) then
554        abort_message='nexca = 0 dans interfoce_cpl'
555        call abort_gcm(modname,abort_message,1)
556      endif
557
558      cumul = .false.
559
560      iloc = maxloc(fder(1:klon))
561      if (check) then
562        if (fder(iloc(1))> 0.) then
563          WRITE(*,*)'**** Debug fder ****'
564          WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))
565        endif
566      endif
567!!$
568!!$      where(fder.gt.0.)
569!!$        fder = 0.
570!!$      endwhere
571
572      call interfoce(itime, dtime, cumul, &
573      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
574      & ocean, npas, nexca, debut, lafin, &
575      & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
576      & fluxlat, fluxsens, fder, albedo, taux, tauy, &
577! -- LOOP
578      & windsp, &
579! -- LOOP
580      & zmasq, &
581      & tsurf_new, alb_new, &
582      & pctsrf_new)
583
584!IM: "slab" ocean
585     else if (ocean == 'slab  ') then
586      DO i=1, knon
587       tsurf_new(i) = tmp_tslab(knindex(i))
588      ENDDO
589      pctsrf_new = tmp_pctsrf_slab
590!
591    else                              ! lecture conditions limites
592      call interfoce(itime, dtime, jour, &
593     &  klon, nisurf, knon, knindex, &
594     &  debut, &
595     &  tsurf_new, pctsrf_new)
596
597    endif
598
599    tsurf_temp = tsurf_new
600    cal = 0.
601    beta = 1.
602    dif_grnd = 0.
603    alb_neig(:) = 0.
604    agesno(:) = 0.
605
606    call calcul_fluxs( klon, knon, nisurf, dtime, &
607     &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
608     &   precip_rain, precip_snow, snow, qsurf,  &
609     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
610     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
611     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
612
613    fder_prev = fder   
614    fder = fder_prev + dflux_s + dflux_l
615
616      iloc = maxloc(fder(1:klon))
617        if (check.and.fder(iloc(1))> 0.) then
618          WRITE(*,*)'**** Debug fder****'
619          WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))
620          WRITE(*,*)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &
621     &                        dflux_s(iloc(1)), dflux_l(iloc(1))
622        endif
623!!$
624!!$      where(fder.gt.0.)
625!!$        fder = 0.
626!!$      endwhere
627
628!IM: flux ocean-atmosphere utile pour le "slab" ocean
629    DO i=1, knon
630     zx_sl(i) = RLVTT
631     if (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
632!IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
633     flux_o(i) = fluxsens(i) + fluxlat(i)
634     tmp_flux_o(knindex(i)) = flux_o(i)
635     tmp_radsol(knindex(i))=radsol(i)
636    ENDDO
637!
638! 2eme appel a interfoce pour le cumul des champs (en particulier
639! fluxsens et fluxlat calcules dans calcul_fluxs)
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, &
650! -- LOOP
651      & windsp, &
652! -- LOOP
653      & zmasq, &
654      & tsurf_new, alb_new, &
655      & pctsrf_new)
656
657    endif
658
659!
660! calcul albedo
661!
662
663    if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then
664      CALL alboc(FLOAT(jour),rlat,alb_eau)
665    else  ! cycle diurne
666      CALL alboc_cd(rmu0,alb_eau)
667    endif
668    DO ii =1, knon
669      alb_new(ii) = alb_eau(knindex(ii))
670    enddo
671
672    z0_new = sqrt(rugos**2 + rugoro**2)
673    alblw(1:knon) = alb_new(1:knon)
674
675!
676  else if (nisurf == is_sic) then
677
678    if (check) write(*,*)'sea ice, nisurf = ',nisurf,'knon=',knon
679!
680! Surface "glace de mer" appel a l'interface avec l'ocean
681!
682!
683    if (ocean == 'couple') then
684
685      cumul =.false.
686
687      iloc = maxloc(fder(1:klon))
688      if (check.and.fder(iloc(1))> 0.) then
689        WRITE(*,*)'**** Debug fder ****'
690        WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))
691      endif
692!!$
693!!$      where(fder.gt.0.)
694!!$        fder = 0.
695!!$      endwhere
696
697      call interfoce(itime, dtime, cumul, &
698      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
699      & ocean, npas, nexca, debut, lafin, &
700      & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
701      & fluxlat, fluxsens, fder, albedo, taux, tauy, &
702! -- LOOP
703      & windsp, &
704! -- LOOP
705      & zmasq, &
706      & tsurf_new, alb_new, &
707      & pctsrf_new)
708
709      tsurf_temp = tsurf_new
710      cal = 0.
711      dif_grnd = 0.
712      beta = 1.0
713
714!IM: "slab" ocean
715    else if (ocean == 'slab  ') then
716!IM ajout sicOBSERVE BEG
717   IF ( ok_slab_sicOBS) THEN
718!                              ! lecture conditions limites
719      CALL interfoce(itime, dtime, jour, &
720             &  klon, nisurf, knon, knindex, &
721             &  debut, &
722             &  tsurf_new, pctsrf_new)
723!
724    tmp_pctsrf_slab=pctsrf_new
725    print*,'jour lecture pctsrf_new sic =',jour
726!
727   ELSE !ok_slab_sicOBS
728     pctsrf_new=tmp_pctsrf_slab
729   ENDIF
730!IM ajout sicOBSERVE END
731!
732     DO ii = 1, knon
733      tsurf_new(ii) = tsurf(ii)
734      IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then
735       snow(ii) = 0.0
736       tsurf_new(ii) = RTT - 1.8
737       IF (soil_model) tsoil(ii,:) = RTT -1.8
738      ENDIF
739     ENDDO
740
741     CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
742     
743     IF (soil_model) THEN
744      CALL soil(dtime, nisurf, knon,snow, tsurf_new, tsoil,soilcap, soilflux)
745      cal(1:knon) = RCPD / soilcap(1:knon)
746      radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
747     ELSE
748      dif_grnd = 1.0 / tau_gl
749      cal = RCPD * calice
750      WHERE (snow > 0.0) cal = RCPD * calsno
751     ENDIF
752     tsurf_temp = tsurf_new
753     beta = 1.0
754!
755    ELSE
756!                              ! lecture conditions limites
757      CALL interfoce(itime, dtime, jour, &
758             &  klon, nisurf, knon, knindex, &
759             &  debut, &
760             &  tsurf_new, pctsrf_new)
761
762!IM cf LF
763      DO ii = 1, knon
764       tsurf_new(ii) = tsurf(ii)
765       IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then
766          snow(ii) = 0.0
767!IM cf LF/JLD         tsurf(ii) = RTT - 1.8
768          tsurf_new(ii) = RTT - 1.8
769          IF (soil_model) tsoil(ii,:) = RTT -1.8
770        endif
771      enddo
772
773      CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
774     
775      IF (soil_model) THEN
776!IM cf LF/JLD        CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)
777         CALL soil(dtime, nisurf, knon,snow, tsurf_new, tsoil,soilcap, soilflux)
778         cal(1:knon) = RCPD / soilcap(1:knon)
779         radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
780         dif_grnd = 1.0 / tau_gl
781      ELSE
782         dif_grnd = 1.0 / tau_gl
783         cal = RCPD * calice
784         WHERE (snow > 0.0) cal = RCPD * calsno
785      ENDIF
786!IMbadtsurf_temp = tsurf
787      tsurf_temp = tsurf_new
788      beta = 1.0
789    ENDIF !ocean ==
790
791    CALL calcul_fluxs( klon, knon, nisurf, dtime, &
792         &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
793         &   precip_rain, precip_snow, snow, qsurf,  &
794         &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
795         &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
796         &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
797!
798    IF (ocean /= 'couple') THEN
799      CALL fonte_neige( klon, knon, nisurf, dtime, &
800             &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
801             &   precip_rain, precip_snow, snow, qsol,  &
802             &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
803             &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
804             &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
805             &   fqcalving,fqfonte,ffonte, run_off_lic_0)
806
807!     calcul albedo
808
809      CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
810      WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
811      zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
812      alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
813     &                    0.6 * (1.0-zfra(1:knon))
814!!      alb_new(1 : knon) = 0.6
815    ENDIF
816!
817!IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
818!
819    DO i = 1, knon
820!
821!IM: faire dependre le coefficient de conduction de la glace de mer
822!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
823!    actuel correspond a 3m de glace de mer, cf. L.Li
824!
825       IF(1.EQ.0) THEN
826        IF(siceh(i).GT.0.) THEN
827         new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)
828        ELSE
829         new_dif_grnd(i) = 0.
830        ENDIF
831       ENDIF !(1.EQ.0) THEN
832!
833     IF (cal(i).GT.1.0e-15) THEN
834       flux_g(i)=(tsurf_new(i)-t_grnd) &
835    &                          * dif_grnd(i) *RCPD/cal(i)
836!   &                          * new_dif_grnd(i) *RCPD/cal(i)
837     ENDIF
838     tmp_flux_g(knindex(i))=flux_g(i)
839!
840!IM: Attention: ne pas initialiser le tmp_radsol puisque c'est deja fait sur is_oce;
841!IM:            tmp_radsol doit etre le flux solaire qui arrive sur l'ocean
842!IM:            et non pas celui qui arrive sur la glace de mer
843!
844    ENDDO
845
846    fder_prev = fder   
847    fder = fder_prev + dflux_s + dflux_l
848
849      iloc = maxloc(fder(1:klon))
850      if (check.and.fder(iloc(1))> 0.) then
851        WRITE(*,*)'**** Debug fder ****'
852        WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))
853        WRITE(*,*)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &
854     &                        dflux_s(iloc(1)), dflux_l(iloc(1))
855      endif
856!!$      where(fder.gt.0.)
857!!$        fder = 0.
858!!$      endwhere
859
860!
861! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
862!
863    if (ocean == 'couple') then
864
865      cumul =.true.
866
867      call interfoce(itime, dtime, cumul, &
868      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
869      & ocean, npas, nexca, debut, lafin, &
870      & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &
871      & fluxlat, fluxsens, fder, albedo, taux, tauy, &
872! -- LOOP
873      & windsp, &
874! -- LOOP
875      & zmasq, &
876      & tsurf_new, alb_new, &
877      & pctsrf_new)
878
879!IM: "slab" ocean
880    else if (ocean == 'slab  ') then
881!
882   IF (check) THEN
883    amn=MIN(tmp_tslab(1),1000.)
884    amx=MAX(tmp_tslab(1),-1000.)
885    DO i=2, klon
886     amn=MIN(tmp_tslab(i),amn)
887     amx=MAX(tmp_tslab(i),amx)
888    ENDDO
889!
890    PRINT*,' debut avant interfoce_slab min max tmp_tslab',amn,amx
891   ENDIF !(check) THEN
892!
893       cumul = .true.
894       tslab = tmp_tslab   
895       call interfoce(klon, debut, itime, dtime, jour, &
896     & tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
897     & tslab, seaice, pctsrf_new)
898!
899       tmp_seaice=seaice
900       tmp_pctsrf_slab=pctsrf_new
901       DO i=1, knon
902         tmp_tslab(knindex(i))=tslab(knindex(i))
903       ENDDO !i
904!
905
906    endif
907
908     
909    z0_new = 0.002
910    z0_new = SQRT(z0_new**2+rugoro**2)
911    alblw(1:knon) = alb_new(1:knon)
912
913
914  else if (nisurf == is_lic) then
915
916    if (check) write(*,*)'glacier, nisurf = ',nisurf
917
918    if (.not. allocated(run_off_lic)) then
919      allocate(run_off_lic(klon), stat = error)
920      if (error /= 0) then
921        abort_message='Pb allocation run_off_lic'
922        call abort_gcm(modname,abort_message,1)
923      endif
924      run_off_lic = 0.
925    endif
926!
927! Surface "glacier continentaux" appel a l'interface avec le sol
928!
929!    call interfsol(nisurf)
930    IF (soil_model) THEN
931        CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil,soilcap, soilflux)
932        cal(1:knon) = RCPD / soilcap(1:knon)
933        radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
934    ELSE
935        cal = RCPD * calice
936        WHERE (snow > 0.0) cal = RCPD * calsno
937    ENDIF
938    beta = 1.0
939    dif_grnd = 0.0
940
941    call calcul_fluxs( klon, knon, nisurf, dtime, &
942     &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
943     &   precip_rain, precip_snow, snow, qsurf,  &
944     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
945     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
946     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
947
948    call fonte_neige( klon, knon, nisurf, dtime, &
949     &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
950     &   precip_rain, precip_snow, snow, qsol,  &
951     &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
952     &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
953     &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
954     &   fqcalving,fqfonte,ffonte, run_off_lic_0)
955
956! passage du run-off des glaciers calcule dans fonte_neige au coupleur
957    bidule=0.
958    bidule(1:knon)= run_off_lic(1:knon)   
959    call gath2cpl(bidule, tmp_rlic, klon, knon,iim,jjm,knindex)
960!
961! calcul albedo
962!
963     CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
964     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
965     zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
966     alb_new(1 : knon)  = alb_neig(1 : knon)*zfra(1:knon) + &
967    &                     0.6 * (1.0-zfra(1:knon))
968!
969!IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
970!       alb_new(1 : knon)  = 0.6 !IM cf FH/GK
971!       alb_new(1 : knon)  = 0.82
972!       alb_new(1 : knon)  = 0.77 !211003 Ksta0.77
973!       alb_new(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5
974!IM: KstaTER0.77 & LMD_ARMIP6   
975        alb_new(1 : knon)  = 0.77
976
977!
978! Rugosite
979!
980    z0_new = rugoro
981!
982! Remplissage des pourcentages de surface
983!
984    pctsrf_new(:,nisurf) = pctsrf(:,nisurf)
985
986    alblw(1:knon) = alb_new(1:knon)
987  else
988    write(*,*)'Index surface = ',nisurf
989    abort_message = 'Index surface non valable'
990    call abort_gcm(modname,abort_message,1)
991  endif
992
993  END SUBROUTINE interfsurf_hq
994
995!
996!#########################################################################
997!
998  SUBROUTINE interfsurf_vent(nisurf, knon   &         
999  &                     )
1000!
1001! Cette routine sert d'aiguillage entre l'atmosphere et la surface en general
1002! (sols continentaux, oceans, glaces) pour les tensions de vents.
1003! En pratique l'interface se fait entre la couche limite du modele
1004! atmospherique (clmain.F) et les routines de surface (sechiba, oasis, ...)
1005!
1006!
1007! L.Fairhead 02/2000
1008!
1009! input:
1010!   nisurf       index de la surface a traiter (1 = sol continental)
1011!   knon         nombre de points de la surface a traiter
1012
1013! Parametres d'entree
1014  integer, intent(IN) :: nisurf
1015  integer, intent(IN) :: knon
1016
1017
1018  return
1019  END SUBROUTINE interfsurf_vent
1020!
1021!#########################################################################
1022!
1023#ifdef CPP_VEGET
1024  SUBROUTINE interfsol(itime, klon, dtime, date0, nisurf, knon, &
1025     & knindex, rlon, rlat, cufi, cvfi, iim, jjm, pctsrf, &
1026     & debut, lafin, ok_veget, &
1027     & plev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
1028     & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1029     & precip_rain, precip_snow, lwdown, swnet, swdown, &
1030     & tsurf, p1lay, ps, radsol, &
1031     & evap, fluxsens, fluxlat, &             
1032     & tsol_rad, tsurf_new, alb_new, alblw, &
1033     & emis_new, z0_new, dflux_l, dflux_s, qsurf)
1034
1035  USE intersurf
1036  USE dimphy, klon_x=>klon
1037  IMPLICIT NONE
1038! Cette routine sert d'interface entre le modele atmospherique et le
1039! modele de sol continental. Appel a sechiba
1040!
1041! L. Fairhead 02/2000
1042!
1043! input:
1044!   itime        numero du pas de temps
1045!   klon         nombre total de points de grille
1046!   dtime        pas de temps de la physique (en s)
1047!   nisurf       index de la surface a traiter (1 = sol continental)
1048!   knon         nombre de points de la surface a traiter
1049!   knindex      index des points de la surface a traiter
1050!   rlon         longitudes de la grille entiere
1051!   rlat         latitudes de la grille entiere
1052!   pctsrf       tableau des fractions de surface de chaque maille
1053!   debut        logical: 1er appel a la physique (lire les restart)
1054!   lafin        logical: dernier appel a la physique (ecrire les restart)
1055!   ok_veget     logical: appel ou non au schema de surface continental
1056!                     (si false calcul simplifie des fluxs sur les continents)
1057!   plev         hauteur de la premiere couche (Pa)     
1058!   u1_lay       vitesse u 1ere couche
1059!   v1_lay       vitesse v 1ere couche
1060!   temp_air     temperature de l'air 1ere couche
1061!   spechum      humidite specifique 1ere couche
1062!   epot_air     temp pot de l'air
1063!   ccanopy      concentration CO2 canopee
1064!   tq_cdrag     cdrag
1065!   petAcoef     coeff. A de la resolution de la CL pour t
1066!   peqAcoef     coeff. A de la resolution de la CL pour q
1067!   petBcoef     coeff. B de la resolution de la CL pour t
1068!   peqBcoef     coeff. B de la resolution de la CL pour q
1069!   precip_rain  precipitation liquide
1070!   precip_snow  precipitation solide
1071!   lwdown       flux IR descendant a la surface
1072!   swnet        flux solaire net
1073!   swdown       flux solaire entrant a la surface
1074!   tsurf        temperature de surface
1075!   p1lay        pression 1er niveau (milieu de couche)
1076!   ps           pression au sol
1077!   radsol       rayonnement net aus sol (LW + SW)
1078!   
1079!
1080! input/output
1081!   run_off      ruissellement total
1082!
1083! output:
1084!   evap         evaporation totale
1085!   fluxsens     flux de chaleur sensible
1086!   fluxlat      flux de chaleur latente
1087!   tsol_rad     
1088!   tsurf_new    temperature au sol
1089!   alb_new      albedo
1090!   emis_new     emissivite
1091!   z0_new       surface roughness
1092!   qsurf        air moisture at surface
1093
1094! Parametres d'entree
1095  integer, intent(IN) :: itime
1096  integer, intent(IN) :: klon
1097  real, intent(IN)    :: dtime
1098  real, intent(IN)    :: date0
1099  integer, intent(IN) :: nisurf
1100  integer, intent(IN) :: knon
1101  integer, intent(IN) :: iim, jjm
1102  integer, dimension(klon), intent(IN) :: knindex
1103  logical, intent(IN) :: debut, lafin, ok_veget
1104  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
1105  real, dimension(klon), intent(IN) :: rlon, rlat
1106  real, dimension(klon), intent(IN) :: cufi, cvfi
1107  real, dimension(klon), intent(IN) :: plev
1108  real, dimension(klon), intent(IN) :: u1_lay, v1_lay
1109  real, dimension(klon), intent(IN) :: temp_air, spechum
1110  real, dimension(klon), intent(IN) :: epot_air, ccanopy
1111  real, dimension(klon), intent(INOUT) :: tq_cdrag
1112  real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
1113  real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
1114  real, dimension(klon), intent(IN) :: precip_rain, precip_snow
1115  real, dimension(klon), intent(IN) :: lwdown, swnet, swdown, ps
1116!IM cf. JP +++
1117  real, dimension(klon) :: swdown_vrai
1118!IM cf. JP ---
1119  real, dimension(klon), intent(IN) :: tsurf, p1lay
1120  real, dimension(klon), intent(IN) :: radsol
1121! Parametres de sortie
1122  real, dimension(klon), intent(OUT):: evap, fluxsens, fluxlat, qsurf
1123  real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new, alblw
1124  real, dimension(klon), intent(OUT):: emis_new, z0_new
1125  real, dimension(klon), intent(OUT):: dflux_s, dflux_l
1126
1127! Local
1128!
1129  integer              :: ii, ij, jj, igrid, ireal, i, index, iglob
1130  integer              :: error
1131  character (len = 20) :: modname = 'interfsol'
1132  character (len = 80) :: abort_message
1133  logical,save              :: check = .FALSE.
1134  real, dimension(klon) :: cal, beta, dif_grnd, capsol
1135! type de couplage dans sechiba
1136!  character (len=10)   :: coupling = 'implicit'
1137! drapeaux controlant les appels dans SECHIBA
1138!  type(control_type), save   :: control_in
1139! Preserved albedo
1140!IM cf. JP +++
1141  real, allocatable, dimension(:), save :: albedo_keep, zlev
1142!IM cf. JP ---
1143! coordonnees geographiques
1144  real, allocatable, dimension(:,:), save :: lalo
1145! pts voisins
1146  integer,allocatable, dimension(:,:), save :: neighbours
1147! fractions continents
1148  real,allocatable, dimension(:), save :: contfrac
1149! resolution de la grille
1150  real, allocatable, dimension (:,:), save :: resolution
1151! correspondance point n -> indices (i,j)
1152  integer, allocatable, dimension(:,:), save :: correspond
1153! offset pour calculer les point voisins
1154  integer, dimension(8,3), save :: off_ini
1155 ! Identifieurs des fichiers restart et histoire
1156  integer, save          :: rest_id, hist_id
1157  integer, save          :: rest_id_stom, hist_id_stom
1158!
1159  real, allocatable, dimension (:,:), save :: lon_scat, lat_scat 
1160
1161  logical, save          :: lrestart_read = .true. , lrestart_write = .false.
1162
1163  real, dimension(klon):: snow
1164  real, dimension(knon,2) :: albedo_out
1165! Pb de nomenclature
1166  real, dimension(klon) :: petA_orc, peqA_orc
1167  real, dimension(klon) :: petB_orc, peqB_orc
1168! Pb de correspondances de grilles
1169  integer, dimension(:), save, allocatable :: ig, jg
1170  integer :: indi, indj
1171  integer, save, allocatable,dimension(:) :: ktindex
1172  REAL, dimension(klon) :: bidule
1173! Essai cdrag
1174  real, dimension(klon) :: cdrag
1175  integer :: jjb,jje,ijb,ije
1176  INTEGER,SAVE :: offset
1177  REAL, dimension(klon2) :: rlon_g,rlat_g
1178  INTEGER, SAVE          :: orch_comm
1179#include "temps.inc"
1180#include "YOMCST.inc"
1181#include "iniprint.h"
1182
1183  if (check) write(lunout,*)'Entree ', modname
1184  if (check) write(lunout,*)'ok_veget = ',ok_veget
1185
1186
1187 
1188! initialisation
1189 
1190  if (debut) then
1191    ALLOCATE(ktindex(klon))
1192  IF ( .NOT. allocated(albedo_keep)) THEN
1193     ALLOCATE(albedo_keep(klon))
1194     ALLOCATE(zlev(klon))
1195  ENDIF
1196! Pb de correspondances de grilles
1197   allocate(ig(klon))
1198   allocate(jg(klon))
1199   ig(1) = 1
1200   jg(1) = 1
1201   indi = 0
1202   indj = 2
1203   do igrid = 2, klon - 1
1204     indi = indi + 1
1205     if ( indi > iim) then
1206       indi = 1
1207       indj = indj + 1
1208     endif
1209     ig(igrid) = indi
1210     jg(igrid) = indj
1211   enddo
1212   ig(klon) = 1
1213   jg(klon) = jjm + 1
1214
1215    if ((.not. allocated(lalo))) then
1216      allocate(lalo(knon,2), stat = error)
1217      if (error /= 0) then
1218        abort_message='Pb allocation lalo'
1219        call abort_gcm(modname,abort_message,1)
1220      endif     
1221    endif
1222    if ((.not. allocated(lon_scat))) then
1223      allocate(lon_scat(iim,jjm+1), stat = error)
1224      if (error /= 0) then
1225        abort_message='Pb allocation lon_scat'
1226        call abort_gcm(modname,abort_message,1)
1227      endif     
1228    endif
1229    if ((.not. allocated(lat_scat))) then
1230      allocate(lat_scat(iim,jjm+1), stat = error)
1231      if (error /= 0) then
1232        abort_message='Pb allocation lat_scat'
1233        call abort_gcm(modname,abort_message,1)
1234      endif     
1235    endif
1236    lon_scat = 0.
1237    lat_scat = 0.
1238    do igrid = 1, knon
1239      index = knindex(igrid)
1240      lalo(igrid,2) = rlon(index)
1241      lalo(igrid,1) = rlat(index)
1242
1243    enddo
1244   
1245
1246   
1247    Call GatherField(rlon,rlon_g,1)
1248    Call GatherField(rlat,rlat_g,1)
1249
1250    IF (phy_rank==0) THEN
1251      index = 1
1252      do jj = 2, jjm
1253        do ij = 1, iim
1254          index = index + 1
1255          lon_scat(ij,jj) = rlon_g(index)
1256          lat_scat(ij,jj) = rlat_g(index)
1257        enddo
1258      enddo
1259     lon_scat(:,1) = lon_scat(:,2)
1260     lat_scat(:,1) = rlat_g(1)
1261     lon_scat(:,jjm+1) = lon_scat(:,2)
1262     lat_scat(:,jjm+1) = rlat_g(klon2)
1263   ENDIF
1264   
1265
1266!
1267! Allouer et initialiser le tableau des voisins et des fraction de continents
1268!
1269    if ( (.not.allocated(neighbours))) THEN
1270      allocate(neighbours(knon,8), stat = error)
1271      if (error /= 0) then
1272        abort_message='Pb allocation neighbours'
1273        call abort_gcm(modname,abort_message,1)
1274      endif
1275    endif
1276    neighbours = -1.
1277    if (( .not. allocated(contfrac))) then
1278      allocate(contfrac(knon), stat = error)
1279      if (error /= 0) then
1280        abort_message='Pb allocation contfrac'
1281        call abort_gcm(modname,abort_message,1)
1282      endif     
1283    endif
1284
1285    do igrid = 1, knon
1286      ireal = knindex(igrid)
1287      contfrac(igrid) = pctsrf(ireal,is_ter)
1288    enddo
1289
1290
1291   CALL Init_neighbours(iim,jjm,knon,neighbours,knindex,pctsrf(:,is_ter))
1292
1293!
1294!  Allocation et calcul resolutions
1295    IF ( (.NOT.ALLOCATED(resolution))) THEN
1296      ALLOCATE(resolution(knon,2), stat = error)
1297      if (error /= 0) then
1298        abort_message='Pb allocation resolution'
1299        call abort_gcm(modname,abort_message,1)
1300      endif
1301    ENDIF
1302    do igrid = 1, knon
1303      ij = knindex(igrid)
1304      resolution(igrid,1) = cufi(ij)
1305      resolution(igrid,2) = cvfi(ij)
1306    enddo 
1307
1308  endif                          ! (fin debut)
1309
1310!
1311! Appel a la routine sols continentaux
1312!
1313  if (lafin) lrestart_write = .true.
1314  if (check) write(lunout,*)'lafin ',lafin,lrestart_write
1315
1316  petA_orc = petBcoef * dtime
1317  petB_orc = petAcoef
1318  peqA_orc = peqBcoef * dtime
1319  peqB_orc = peqAcoef
1320
1321  cdrag = 0.
1322  cdrag(1:knon) = tq_cdrag(1:knon)
1323
1324!IM cf. JP +++
1325! zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/287.05*temp_air(1:knon))*9.80665)
1326  zlev(1:knon) = (100.*plev(1:knon))/((ps(1:knon)/RD*temp_air(1:knon))*RG)
1327!IM cf. JP ---
1328
1329
1330! PF et PASB
1331!   where(cdrag > 0.01)
1332!     cdrag = 0.01
1333!   endwhere
1334!  write(*,*)'Cdrag = ',minval(cdrag),maxval(cdrag)
1335
1336!
1337! Init Orchidee
1338!
1339!  if (pole_nord) then
1340!    offset=0
1341!    ktindex(:)=ktindex(:)+iim-1
1342!  else
1343!    offset = klon_begin-1+iim-1
1344!    ktindex(:)=ktindex(:)+MOD(offset,iim)
1345!    offset=offset-MOD(offset,iim)
1346!  endif
1347 
1348  PRINT *,'ORCHIDEE ------> KNON : ',knon
1349 
1350   
1351  if (debut) then
1352    CALL Get_orchidee_communicator(knon,orch_comm)
1353    IF (knon /=0) THEN
1354      CALL Init_orchidee_index(iim,knon,orch_comm,knindex,offset,ktindex)
1355   
1356      call intersurf_main (itime+itau_phy-1, iim, jjm+1,offset, knon, ktindex, &
1357       & orch_comm,dtime, lrestart_read, lrestart_write, lalo, &
1358       & contfrac, neighbours, resolution, date0, &
1359       & zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
1360       & cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
1361       & precip_rain, precip_snow, lwdown, swnet, swdown, ps, &
1362       & evap, fluxsens, fluxlat, coastalflow, riverflow, &
1363       & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
1364       & lon_scat, lat_scat)
1365
1366     ENDIF
1367!IM cf. JP +++
1368    albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
1369!IM cf. JP ---
1370
1371  endif
1372
1373!IM cf. JP +++
1374!  swdown_vrai(1:knon) = swnet(1:knon)/(1. - albedo_keep(1:knon))
1375  swdown_vrai(1:knon) = swdown(1:knon)
1376
1377!IM cf. JP ---
1378    IF (knon /=0) THEN
1379   
1380      call intersurf_main (itime+itau_phy, iim, jjm+1,offset, knon, ktindex, &
1381       & orch_comm,dtime, lrestart_read, lrestart_write, lalo, &
1382       & contfrac, neighbours, resolution, date0, &
1383       & zlev,  u1_lay, v1_lay, spechum, temp_air, epot_air, ccanopy, &
1384       & cdrag, petA_orc, peqA_orc, petB_orc, peqB_orc, &
1385       & precip_rain, precip_snow, lwdown, swnet, swdown_vrai, ps, &
1386       & evap, fluxsens, fluxlat, coastalflow, riverflow, &
1387       & tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
1388       & lon_scat, lat_scat)
1389
1390    ENDIF
1391!IM cf. JP +++
1392    albedo_keep(1:knon) = (albedo_out(1:knon,1)+albedo_out(1:knon,2))/2.
1393!IM cf. JP ---
1394
1395    bidule=0.
1396    bidule(1:knon)=riverflow(1:knon)
1397    call gath2cpl(bidule, tmp_rriv, klon, knon,iim,jjphy_nb,knindex)
1398    bidule=0.
1399    bidule(1:knon)=coastalflow(1:knon)
1400    call gath2cpl(bidule, tmp_rcoa, klon, knon,iim,jjphy_nb,knindex)
1401    alb_new(1:knon) = albedo_out(1:knon,1)
1402    alblw(1:knon) = albedo_out(1:knon,2)
1403
1404
1405! Convention orchidee: positif vers le haut
1406  fluxsens(1:knon) = -1. * fluxsens(1:knon)
1407  fluxlat(1:knon)  = -1. * fluxlat(1:knon)
1408
1409!  evap     = -1. * evap
1410
1411  if (debut) lrestart_read = .false.
1412
1413  END SUBROUTINE interfsol
1414 
1415  SUBROUTINE Init_orchidee_index(iim,knon,orch_comm,knindex,offset,ktindex)
1416  USE dimphy
1417  IMPLICIT NONE
1418    INTEGER,INTENT(IN)  :: iim
1419    INTEGER,INTENT(IN)  :: knon
1420    INTEGER,INTENT(IN)  :: orch_comm
1421    INTEGER,INTENT(IN)  :: knindex(knon)
1422    INTEGER,INTENT(OUT) :: offset
1423    INTEGER,INTENT(OUT) :: ktindex(knon)
1424
1425#ifdef CPP_PARA
1426    INCLUDE 'mpif.h'
1427    INTEGER :: status(MPI_STATUS_SIZE)
1428#endif
1429    INTEGER :: MyLastPoint
1430    INTEGER :: LastPoint
1431    INTEGER :: mpi_rank
1432    INTEGER :: mpi_size
1433    INTEGER :: ierr   
1434   
1435    MyLastPoint=klon_begin-1+knindex(knon)+iim-1
1436
1437    IF (.NOT. monocpu) THEN
1438#ifdef CPP_PARA   
1439      call MPI_COMM_SIZE(orch_comm,mpi_size,ierr)
1440      call MPI_COMM_RANK(orch_comm,mpi_rank,ierr)
1441#endif
1442    ELSE
1443      mpi_rank=0
1444      mpi_size=1
1445    ENDIF
1446   
1447    IF (.NOT. monocpu) THEN
1448      IF (mpi_rank /= 0) then
1449#ifdef CPP_PARA
1450        CALL MPI_RECV(LastPoint,1,MPI_INTEGER,mpi_rank-1,1234,orch_comm,status,ierr)
1451#endif
1452      ENDIF
1453       
1454      IF (mpi_rank /= mpi_size-1) THEN
1455#ifdef CPP_PARA
1456        CALL MPI_SEND(MyLastPoint,1,MPI_INTEGER,mpi_rank+1,1234,orch_comm,ierr) 
1457#endif
1458      ENDIF
1459    ENDIF
1460   
1461    IF (mpi_rank==0) THEN
1462      offset=0
1463    ELSE
1464     offset=LastPoint-MOD(LastPoint,iim)
1465    ENDIF
1466     
1467    ktindex(:)=knindex(:)+(klon_begin+iim-1)-offset-1   
1468   
1469
1470   END SUBROUTINE  Init_orchidee_index
1471
1472 
1473  SUBROUTINE Get_orchidee_communicator(knon,orch_comm)
1474  USE dimphy, only : phy_rank
1475#ifdef CPP_PARA
1476  USE parallel, only : COMM_LMDZ
1477#endif
1478  IMPLICIT NONE
1479#ifdef CPP_PARA
1480    include 'mpif.h'
1481#endif   
1482    INTEGER,INTENT(IN)  :: knon
1483    INTEGER,INTENT(OUT) :: orch_comm
1484   
1485    INTEGER :: color
1486    INTEGER :: ierr
1487   
1488    IF (knon==0) THEN
1489      color = 0
1490    ELSE
1491      color = 1
1492    ENDIF
1493
1494#ifdef CPP_PARA   
1495    CALL MPI_COMM_SPLIT(COMM_LMDZ,color,phy_rank,orch_comm,ierr)
1496#endif
1497   
1498  END SUBROUTINE Get_orchidee_communicator
1499   
1500   
1501  SUBROUTINE Init_neighbours(iim,jjm,knon,neighbours,ktindex,pctsrf)
1502#ifdef CPP_PARA
1503  USE parallel,only : COMM_LMDZ
1504#endif
1505  USE dimphy
1506  IMPLICIT NONE
1507#ifdef CPP_PARA
1508  include 'mpif.h'
1509#endif
1510  INTEGER :: iim,jjm
1511  INTEGER :: knon
1512  INTEGER :: neighbours(knon,8)
1513  INTEGER :: ktindex(knon)
1514  REAL :: pctsrf(klon)
1515 
1516  INTEGER :: knon_nb(0:phy_size-1)
1517  INTEGER,DIMENSION(0:phy_size-1) :: displs,sendcount
1518  INTEGER,ALLOCATABLE :: ktindex_g(:)
1519  REAL*8  :: pctsrf_g(klon2)
1520  INTEGER,ALLOCATABLE ::neighbours_g(:,:)
1521  INTEGER :: knon_g
1522  REAL*8 :: correspond(iim,jjm+1)
1523  INTEGER :: i,igrid,jj,ij,iglob,ierr,ireal,index
1524  integer, dimension(8,3) :: off_ini
1525  integer, dimension(8)   :: offset 
1526  INTEGER :: ktindex_p(knon)
1527
1528  IF (monocpu) THEN
1529    knon_nb(:)=knon
1530  ELSE 
1531
1532#ifdef CPP_PARA 
1533    CALL MPI_GATHER(knon,1,MPI_INTEGER,knon_nb,1,MPI_INTEGER,0,COMM_LMDZ,ierr)
1534#endif
1535 
1536  ENDIF
1537   
1538  IF (phy_rank==0) THEN
1539    knon_g=sum(knon_nb(:))
1540    ALLOCATE(ktindex_g(knon_g))
1541    ALLOCATE(neighbours_g(knon_g,8))
1542    neighbours_g(:,:)=-1
1543    displs(0)=0
1544    DO i=1,phy_size-1
1545      displs(i)=displs(i-1)+knon_nb(i-1)
1546    ENDDO 
1547  ENDIF
1548 
1549  ktindex_p(:)=ktindex(:)+klon_begin-1+iim-1
1550
1551  IF (monocpu) THEN
1552    ktindex_g(:)=ktindex_p(:)
1553  ELSE
1554
1555#ifdef CPP_PARA 
1556    CALL MPI_GATHERV(ktindex_p,knon,MPI_INTEGER,ktindex_g,knon_nb,displs,MPI_INTEGER,0,COMM_LMDZ,ierr)
1557#endif
1558
1559  ENDIF
1560     
1561  CALL GatherField(pctsrf,pctsrf_g,1)
1562 
1563  IF (phy_rank==0) THEN
1564!  Initialisation des offset   
1565!
1566! offset bord ouest
1567   off_ini(1,1) = - iim  ; off_ini(2,1) = - iim + 1; off_ini(3,1) = 1
1568   off_ini(4,1) = iim + 1; off_ini(5,1) = iim      ; off_ini(6,1) = 2 * iim - 1
1569   off_ini(7,1) = iim -1 ; off_ini(8,1) = - 1
1570! offset point normal
1571   off_ini(1,2) = - iim  ; off_ini(2,2) = - iim + 1; off_ini(3,2) = 1
1572   off_ini(4,2) = iim + 1; off_ini(5,2) = iim      ; off_ini(6,2) = iim - 1
1573   off_ini(7,2) = -1     ; off_ini(8,2) = - iim - 1
1574! offset bord   est
1575   off_ini(1,3) = - iim; off_ini(2,3) = - 2 * iim + 1; off_ini(3,3) = - iim + 1
1576   off_ini(4,3) =  1   ; off_ini(5,3) = iim          ; off_ini(6,3) = iim - 1
1577   off_ini(7,3) = -1   ; off_ini(8,3) = - iim - 1
1578!
1579!
1580! Attention aux poles
1581!
1582    do igrid = 1, knon_g
1583      index = ktindex_g(igrid)
1584          jj = int((index - 1)/iim) + 1
1585          ij = index - (jj - 1) * iim
1586      correspond(ij,jj) = igrid
1587    enddo
1588
1589    do igrid = 1, knon_g
1590      iglob = ktindex_g(igrid)
1591      if (mod(iglob, iim) == 1) then
1592        offset = off_ini(:,1)
1593      else if(mod(iglob, iim) == 0) then
1594        offset = off_ini(:,3)
1595      else
1596        offset = off_ini(:,2)
1597      endif
1598      do i = 1, 8
1599        index = iglob + offset(i)
1600        ireal = (min(max(1, index - iim + 1), klon2))
1601        if (pctsrf_g(ireal) > EPSFRA) then
1602          jj = int((index - 1)/iim) + 1
1603          ij = index - (jj - 1) * iim
1604            neighbours_g(igrid, i) = correspond(ij, jj)
1605        endif
1606      enddo
1607    enddo
1608
1609!    DO i=0,phy_size-1
1610!      displs(i)=displs(i)*8
1611!      sendcount(i)=knon_nb(i)*8
1612!    ENDDO
1613 
1614  ENDIF
1615 
1616  DO i=1,8
1617    IF (monocpu) THEN
1618      neighbours(:,i)=neighbours_g(:,i)
1619    ELSE
1620#ifdef CPP_PARA
1621    CALL MPI_SCATTERV(neighbours_g(:,i),knon_nb,displs,MPI_INTEGER,neighbours(:,i),knon,MPI_INTEGER,0,COMM_LMDZ,ierr)
1622#endif
1623    ENDIF
1624  ENDDO
1625 
1626  END SUBROUTINE Init_neighbours
1627#endif
1628!
1629!#########################################################################
1630!
1631  SUBROUTINE interfoce_cpl(itime, dtime, cumul, &
1632      & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &
1633      & ocean, npas, nexca, debut, lafin, &
1634      & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &
1635      & fluxlat, fluxsens, fder, albsol, taux, tauy, &
1636! -- LOOP
1637      & windsp, &
1638! -- LOOP
1639      & zmasq, &
1640      & tsurf_new, alb_new, &
1641      & pctsrf_new)
1642
1643   USE ioipsl
1644   USE dimphy, only : jjphy_nb, iiphy_begin,iiphy_end,phy_rank,phy_size, monocpu
1645   USE iophy
1646#ifdef CPP_PARA
1647   USE parallel, only: pole_nord,pole_sud,COMM_LMDZ
1648#endif
1649#ifdef CPP_PSMILE
1650   USE oasis
1651#endif
1652   USE write_field_phy
1653   implicit none
1654#include "indicesol.inc"
1655#include "YOMCST.inc"
1656! Cette routine sert d'interface entre le modele atmospherique et un
1657! coupleur avec un modele d'ocean 'complet' derriere
1658!
1659! Le modele de glace qu'il est prevu d'utiliser etant couple directement a
1660! l'ocean presentement, on va passer deux fois dans cette routine par pas de
1661! temps physique, une fois avec les points oceans et l'autre avec les points
1662! glace. A chaque pas de temps de couplage, la lecture des champs provenant
1663! du coupleur se fera "dans" l'ocean et l'ecriture des champs a envoyer
1664! au coupleur "dans" la glace. Il faut donc des tableaux de travail "tampons"
1665! dimensionnes sur toute la grille qui remplissent les champs sur les
1666! domaines ocean/glace quand il le faut. Il est aussi necessaire que l'index
1667! ocean soit traiter avant l'index glace (sinon tout intervertir)
1668!
1669!
1670! L. Fairhead 02/2000
1671!
1672! input:
1673!   itime        numero du pas de temps
1674!   iim, jjm     nbres de pts de grille
1675!   dtime        pas de temps de la physique
1676!   klon         nombre total de points de grille
1677!   nisurf       index de la surface a traiter (1 = sol continental)
1678!   pctsrf       tableau des fractions de surface de chaque maille
1679!   knon         nombre de points de la surface a traiter
1680!   knindex      index des points de la surface a traiter
1681!   rlon         longitudes
1682!   rlat         latitudes
1683!   debut        logical: 1er appel a la physique
1684!   lafin        logical: dernier appel a la physique
1685!   ocean        type d'ocean
1686!   nexca        frequence de couplage
1687!   swdown       flux solaire entrant a la surface
1688!   lwdown       flux IR net a la surface
1689!   precip_rain  precipitation liquide
1690!   precip_snow  precipitation solide
1691!   evap         evaporation
1692!   tsurf        temperature de surface
1693!   fder         derivee dF/dT
1694!   albsol       albedo du sol (coherent avec swdown)
1695!   taux         tension de vent en x
1696!   tauy         tension de vent en y
1697! -- LOOP
1698!    windsp       module du vent a 10m
1699! -- LOOP
1700!   nexca        frequence de couplage
1701!   zmasq        masque terre/ocean
1702!
1703!
1704! output:
1705!   tsurf_new    temperature au sol
1706!   alb_new      albedo
1707!   pctsrf_new   nouvelle repartition des surfaces
1708!   alb_ice      albedo de la glace
1709!
1710
1711
1712! Parametres d'entree
1713  integer, intent(IN) :: itime
1714  integer, intent(IN) :: iim, jjm
1715  real, intent(IN) :: dtime
1716  integer, intent(IN) :: klon
1717  integer, intent(IN) :: nisurf
1718  integer, intent(IN) :: knon
1719  real, dimension(klon,nbsrf), intent(IN) :: pctsrf
1720  integer, dimension(klon), intent(in) :: knindex
1721  logical, intent(IN) :: debut, lafin
1722  real, dimension(klon), intent(IN) :: rlon, rlat
1723  character (len = 6)  :: ocean
1724  real, dimension(klon), intent(IN) :: lwdown, swdown
1725  real, dimension(klon), intent(IN) :: precip_rain, precip_snow
1726  real, dimension(klon), intent(IN) :: tsurf, fder, albsol, taux, tauy
1727! -- LOOP
1728   real, dimension(klon), intent(IN) :: windsp
1729! -- LOOP
1730  INTEGER              :: nexca, npas, kstep
1731  real, dimension(klon), intent(IN) :: zmasq
1732  real, dimension(klon), intent(IN) :: fluxlat, fluxsens
1733  logical, intent(IN)               :: cumul
1734  real, dimension(klon), intent(INOUT) :: evap
1735
1736! Parametres de sortie
1737  real, dimension(klon), intent(OUT):: tsurf_new, alb_new
1738  real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new
1739
1740! Variables locales
1741  integer                    :: j, error, sum_error, ig, cpl_index,i
1742! -- LOOP
1743  INTEGER :: nsrf
1744! -- LOOP
1745  character (len = 20) :: modname = 'interfoce_cpl'
1746  character (len = 80) :: abort_message
1747  logical,save              :: check = .FALSE.
1748! variables pour moyenner les variables de couplage
1749  real, allocatable, dimension(:,:),save :: cpl_sols, cpl_nsol, cpl_rain
1750  real, allocatable, dimension(:,:),save :: cpl_snow, cpl_evap, cpl_tsol
1751  real, allocatable, dimension(:,:),save :: cpl_fder, cpl_albe, cpl_taux
1752! -- LOOP
1753  real, allocatable, dimension(:,:),save :: cpl_windsp
1754! -- LOOP
1755  real, allocatable, dimension(:,:),save :: cpl_tauy
1756  REAL, ALLOCATABLE, DIMENSION(:,:),SAVE :: cpl_rriv, cpl_rcoa, cpl_rlic
1757!!$
1758! variables tampons avant le passage au coupleur
1759  real, allocatable, dimension(:,:,:),save :: tmp_sols, tmp_nsol, tmp_rain
1760  real, allocatable, dimension(:,:,:),save :: tmp_snow, tmp_evap, tmp_tsol
1761  real, allocatable, dimension(:,:,:),save :: tmp_fder, tmp_albe, tmp_taux
1762! -- LOOP
1763  real, allocatable, dimension(:,:,:),save :: tmp_windsp
1764! -- LOOP
1765!!$  real, allocatable, dimension(:,:,:),save :: tmp_tauy, tmp_rriv, tmp_rcoa
1766  REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: tmp_tauy
1767! variables a passer au coupleur
1768!ym  real, dimension(iim, jjm+1) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice
1769!ym  real, dimension(iim, jjm+1) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
1770!ym  REAL, DIMENSION(iim, jjm+1) :: wri_evap_sea, wri_rcoa, wri_rriv
1771!ym  REAL, DIMENSION(iim, jjm+1) :: wri_rain, wri_snow, wri_taux, wri_tauy
1772!ym  REAL, DIMENSION(iim, jjm+1) :: wri_calv
1773!ym  REAL, DIMENSION(iim, jjm+1) :: wri_tauxx, wri_tauyy, wri_tauzz
1774!ym  REAL, DIMENSION(iim, jjm+1) :: tmp_lon, tmp_lat
1775
1776  real, dimension(iim, jjphy_nb) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice
1777  real, dimension(iim, jjphy_nb) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice
1778  REAL, DIMENSION(iim, jjphy_nb) :: wri_evap_sea, wri_rcoa, wri_rriv
1779  REAL, DIMENSION(iim, jjphy_nb) :: wri_rain, wri_snow, wri_taux, wri_tauy
1780  REAL, DIMENSION(iim, jjphy_nb) :: wri_calv
1781  REAL, DIMENSION(iim, jjphy_nb) :: wri_tauxx, wri_tauyy, wri_tauzz
1782  REAL, DIMENSION(iim, jjphy_nb) :: tmp_lon, tmp_lat
1783  REAL, DIMENSION(iim, jjphy_nb) :: wri_windsp
1784
1785! variables relues par le coupleur
1786! read_sic = fraction de glace
1787! read_sit = temperature de glace
1788  real, allocatable, dimension(:,:),save :: read_sst, read_sic, read_sit
1789  real, allocatable, dimension(:,:),save :: read_alb_sic
1790! variable tampon
1791  real, dimension(klon)       :: tamp_sic
1792! sauvegarde des fractions de surface d'un pas de temps a l'autre apres
1793! l'avoir lu
1794  real, allocatable,dimension(:,:),save :: pctsrf_sav
1795  real, dimension(iim, jjphy_nb, 3) :: tamp_srf
1796  integer, allocatable, dimension(:), save :: tamp_ind
1797  real, allocatable, dimension(:,:),save :: tamp_zmasq
1798  real, dimension(iim, jjphy_nb) :: deno
1799  integer                     :: idtime
1800  integer, allocatable,dimension(:),save :: unity
1801!
1802  logical, save    :: first_appel = .true.
1803  logical,save          :: print
1804!maf
1805! variables pour avoir une sortie IOIPSL des champs echanges
1806  CHARACTER*80,SAVE :: clintocplnam, clfromcplnam
1807  INTEGER, SAVE :: jf,nhoridct,nidct
1808  INTEGER, SAVE :: nhoridcs,nidcs
1809  INTEGER :: ndexct(iim*(jjm+1)),ndexcs(iim*(jjm+1))
1810  REAL :: zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian
1811  INTEGER,save :: idayref
1812!med  integer :: itau_w
1813  integer,save :: itau_w
1814! -- LOOP
1815  integer :: nb_interf_cpl
1816! -- LOOP
1817
1818  real :: Up,Down
1819  integer :: ierr
1820  integer :: il_time_secs
1821  real :: tmp_field(klon)
1822 
1823#include "param_cou.h"
1824#include "inc_cpl.h"
1825#include "temps.inc"
1826#include "iniprint.h"
1827
1828#ifdef CPP_PARA
1829  include 'mpif.h'
1830  integer :: status(MPI_STATUS_SIZE)
1831#endif
1832
1833!
1834! Initialisation
1835!
1836  if (check) write(*,*)'Entree ',modname,'nisurf = ',nisurf
1837 
1838  if (first_appel) then
1839    error = 0
1840    allocate(unity(klon), stat = error)
1841    if ( error  /=0) then
1842      abort_message='Pb allocation variable unity'
1843      call abort_gcm(modname,abort_message,1)
1844    endif
1845    allocate(pctsrf_sav(klon,nbsrf), stat = error)
1846    if ( error  /=0) then
1847      abort_message='Pb allocation variable pctsrf_sav'
1848      call abort_gcm(modname,abort_message,1)
1849    endif
1850    pctsrf_sav = 0.
1851
1852    do ig = 1, klon
1853      unity(ig) = ig
1854    enddo
1855    sum_error = 0
1856    allocate(cpl_sols(klon,2), stat = error); sum_error = sum_error + error
1857    allocate(cpl_nsol(klon,2), stat = error); sum_error = sum_error + error
1858    allocate(cpl_rain(klon,2), stat = error); sum_error = sum_error + error
1859    allocate(cpl_snow(klon,2), stat = error); sum_error = sum_error + error
1860    allocate(cpl_evap(klon,2), stat = error); sum_error = sum_error + error
1861    allocate(cpl_tsol(klon,2), stat = error); sum_error = sum_error + error
1862    allocate(cpl_fder(klon,2), stat = error); sum_error = sum_error + error
1863    allocate(cpl_albe(klon,2), stat = error); sum_error = sum_error + error
1864    allocate(cpl_taux(klon,2), stat = error); sum_error = sum_error + error
1865! -- LOOP
1866     allocate(cpl_windsp(klon,2), stat = error); sum_error = sum_error + error
1867! -- LOOP
1868    allocate(cpl_tauy(klon,2), stat = error); sum_error = sum_error + error
1869!ym    ALLOCATE(cpl_rriv(iim,jjm+1), stat=error); sum_error = sum_error + error
1870!ym    ALLOCATE(cpl_rcoa(iim,jjm+1), stat=error); sum_error = sum_error + error
1871!ym    ALLOCATE(cpl_rlic(iim,jjm+1), stat=error); sum_error = sum_error + error
1872    ALLOCATE(cpl_rriv(iim,jjphy_nb), stat=error); sum_error = sum_error + error
1873    ALLOCATE(cpl_rcoa(iim,jjphy_nb), stat=error); sum_error = sum_error + error
1874    ALLOCATE(cpl_rlic(iim,jjphy_nb), stat=error); sum_error = sum_error + error
1875
1876
1877!!
1878!ym    allocate(read_sst(iim, jjm+1), stat = error); sum_error = sum_error + error
1879!ym    allocate(read_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
1880!ym    allocate(read_sit(iim, jjm+1), stat = error); sum_error = sum_error + error
1881!ym    allocate(read_alb_sic(iim, jjm+1), stat = error); sum_error = sum_error + error
1882    allocate(read_sst(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1883    allocate(read_sic(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1884    allocate(read_sit(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1885    allocate(read_alb_sic(iim, jjphy_nb), stat = error); sum_error = sum_error + error
1886    read_sst=0.
1887    read_sic=0.
1888    read_sit=0.
1889    read_alb_sic=0.
1890    if (sum_error /= 0) then
1891      abort_message='Pb allocation variables couplees'
1892      call abort_gcm(modname,abort_message,1)
1893    endif
1894    cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
1895    cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
1896    cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.; cpl_rlic = 0.
1897! -- LOOP
1898     cpl_windsp = 0.
1899! -- LOOP
1900
1901    sum_error = 0
1902    allocate(tamp_ind(klon), stat = error); sum_error = sum_error + error
1903!ym    allocate(tamp_zmasq(iim, jjm+1), stat = error); sum_error = sum_error + error   
1904    allocate(tamp_zmasq(iim, jjphy_nb), stat = error); sum_error = sum_error + error   
1905    tamp_zmasq=1.
1906   
1907    do ig = 1, klon
1908      tamp_ind(ig) = ig
1909    enddo
1910    call gath2cpl(zmasq, tamp_zmasq, klon, klon, iim, jjphy_nb, tamp_ind)
1911!
1912! initialisation couplage
1913!
1914    idtime = int(dtime)
1915#ifdef CPP_COUPLE
1916#ifdef CPP_PSMILE
1917   CALL inicma(iim, (jjm+1))
1918#else
1919   if (.not. monocpu) then
1920      abort_message='coupleur parallele uniquement avec PSMILE'
1921      call abort_gcm(modname,abort_message,1)
1922   endif
1923   call inicma(npas , nexca, idtime,(jjm+1)*iim)
1924#endif
1925#endif
1926!
1927! initialisation sorties netcdf
1928!
1929 !ym  IO de check deconnecte pour le moment en //
1930    IF (monocpu) THEN
1931    idayref = day_ini
1932    CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)
1933    CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)
1934    DO i = 1, iim
1935      zx_lon(i,1) = rlon(i+1)
1936      zx_lon(i,jjm+1) = rlon(i+1)
1937    ENDDO
1938    CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)
1939    clintocplnam="cpl_atm_tauflx"
1940    CALL histbeg(clintocplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
1941       & itau_phy,zjulian,dtime,nhoridct,nidct)
1942! no vertical axis
1943    CALL histdef(nidct, 'tauxe','tauxe', &
1944         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1945    CALL histdef(nidct, 'tauyn','tauyn', &
1946         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1947    CALL histdef(nidct, 'tmp_lon','tmp_lon', &
1948         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1949    CALL histdef(nidct, 'tmp_lat','tmp_lat', &
1950         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1951    DO jf=1,jpflda2o1 + jpflda2o2
1952      CALL histdef(nidct, cl_writ(jf),cl_writ(jf), &
1953         & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1954    END DO
1955    CALL histend(nidct)
1956    CALL histsync(nidct)
1957
1958    clfromcplnam="cpl_atm_sst"
1959    CALL histbeg(clfromcplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &
1960       & 0,zjulian,dtime,nhoridcs,nidcs)
1961! no vertical axis
1962    DO jf=1,jpfldo2a
1963      CALL histdef(nidcs, cl_read(jf),cl_read(jf), &
1964         & "-",iim, jjm+1, nhoridcs, 1, 1, 1, -99, 32, "inst", dtime,dtime)
1965    END DO
1966    CALL histend(nidcs)
1967    CALL histsync(nidcs)
1968
1969    ENDIF    ! monocpu
1970   
1971! pour simuler la fonte des glaciers antarctiques
1972!
1973!ym => pour le moment, c'est en commentaire, donc je squizze
1974
1975!ym    surf_maille = (4. * rpi * ra**2) / (iim * (jjm +1))
1976!ym    ALLOCATE(coeff_iceberg(iim,jjm+1), stat=error)
1977!ym    if (error /= 0) then
1978!ym      abort_message='Pb allocation variable coeff_iceberg'
1979!ym      call abort_gcm(modname,abort_message,1)
1980!ym    endif
1981!ym    open (12,file='flux_iceberg',form='formatted',status='old')
1982!ym    read (12,*) coeff_iceberg
1983!ym    close (12)
1984!ym    num_antarctic = max(1, count(coeff_iceberg > 0))
1985   
1986    first_appel = .false.
1987  endif ! fin if (first_appel)
1988
1989! Initialisations
1990
1991! calcul des fluxs a passer
1992! -- LOOP
1993  nb_interf_cpl = nb_interf_cpl + 1 
1994  if (check) write(lunout,*)'passage dans interface_surf.F90 :  ',nb_interf_cpl
1995! -- LOOP
1996  cpl_index = 1
1997  if (nisurf == is_sic) cpl_index = 2
1998  if (cumul) then
1999! -- LOOP
2000      if (check) write(lunout,*)'passage dans cumul '
2001      if (check) write(lunout,*)'valeur de cpl_index ', cpl_index
2002! -- LOOP 
2003    if (check) write(*,*) modname, 'cumul des champs'
2004    do ig = 1, knon
2005      cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) &
2006       &                          + swdown(ig)      / FLOAT(nexca)
2007      cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) &
2008       &                          + (lwdown(ig) + fluxlat(ig) +fluxsens(ig))&
2009       &                                / FLOAT(nexca)
2010      cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) &
2011       &                          + precip_rain(ig) / FLOAT(nexca)
2012      cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) &
2013       &                          + precip_snow(ig) / FLOAT(nexca)
2014      cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) &
2015       &                          + evap(ig)        / FLOAT(nexca)
2016      cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) &
2017       &                          + tsurf(ig)       / FLOAT(nexca)
2018      cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) &
2019       &                          + fder(ig)        / FLOAT(nexca)
2020      cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) &
2021       &                          + albsol(ig)      / FLOAT(nexca)
2022      cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) &
2023       &                          + taux(ig)        / FLOAT(nexca)
2024      cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) &
2025       &                          + tauy(ig)        / FLOAT(nexca)
2026! -- LOOP
2027      IF (cpl_index .EQ. 1) THEN
2028      cpl_windsp(ig,cpl_index) = cpl_windsp(ig,cpl_index) &
2029       &                          + windsp(ig)      / FLOAT(nexca)
2030      ENDIF
2031! -- LOOP
2032    enddo
2033    IF (cpl_index .EQ. 1) THEN
2034        cpl_rriv(:,:) = cpl_rriv(:,:) + tmp_rriv(:,:) / FLOAT(nexca)
2035        cpl_rcoa(:,:) = cpl_rcoa(:,:) + tmp_rcoa(:,:) / FLOAT(nexca)
2036        cpl_rlic(:,:) = cpl_rlic(:,:) + tmp_rlic(:,:) / FLOAT(nexca)
2037    ENDIF
2038  endif
2039
2040  if (mod(itime, nexca) == 1) then
2041!
2042! Demande des champs au coupleur
2043!
2044! Si le domaine considere est l'ocean, on lit les champs venant du coupleur
2045!
2046    if (nisurf == is_oce .and. .not. cumul) then
2047      if (check) write(*,*)'rentree fromcpl, itime-1 = ',itime-1
2048#ifdef CPP_COUPLE
2049#ifdef CPP_PSMILE
2050      il_time_secs=(itime-1)*dtime
2051      CALL fromcpl(il_time_secs, iim, jjphy_nb,                           &
2052     &        read_sst, read_sic, read_sit, read_alb_sic)
2053#else
2054     if (.not. monocpu) then
2055        abort_message='coupleur parallele uniquement avec PSMILE'
2056        call abort_gcm(modname,abort_message,1)
2057     endif
2058
2059      call fromcpl(itime-1,(jjm+1)*iim,                                  &
2060     &        read_sst, read_sic, read_sit, read_alb_sic)
2061#endif
2062#endif
2063!
2064! sorties NETCDF des champs recus
2065!
2066     if (monocpu) THEN
2067       ndexcs(:)=0
2068       itau_w = itau_phy + itime
2069       CALL histwrite(nidcs,cl_read(1),itau_w,read_sst,iim*(jjm+1),ndexcs)
2070       CALL histwrite(nidcs,cl_read(2),itau_w,read_sic,iim*(jjm+1),ndexcs)
2071       CALL histwrite(nidcs,cl_read(3),itau_w,read_alb_sic,iim*(jjm+1),ndexcs)
2072       CALL histwrite(nidcs,cl_read(4),itau_w,read_sit,iim*(jjm+1),ndexcs)
2073       CALL histsync(nidcs)
2074     endif
2075! pas utile      IF (npas-itime.LT.nexca )CALL histclo(nidcs)
2076
2077!ym      do j = 1, jjm + 1
2078       do j = 1, jjphy_nb
2079         do ig = 1, iim
2080          if (abs(1. - read_sic(ig,j)) < 0.00001) then
2081            read_sst(ig,j) = RTT - 1.8
2082            read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
2083            read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
2084          else if (abs(read_sic(ig,j)) < 0.00001) then
2085            read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
2086            read_sit(ig,j) = read_sst(ig,j)
2087            read_alb_sic(ig,j) =  0.6
2088          else
2089            read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))
2090            read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)
2091            read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)
2092          endif
2093        enddo
2094      enddo
2095!
2096! transformer read_sic en pctsrf_sav
2097!
2098      call cpl2gath(read_sic, tamp_sic , klon, klon,iim,jjphy_nb, unity)
2099      do ig = 1, klon
2100        IF (pctsrf(ig,is_oce) > epsfra .OR.            &
2101     &             pctsrf(ig,is_sic) > epsfra) THEN
2102          pctsrf_sav(ig,is_sic) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
2103     &                               * tamp_sic(ig)
2104          pctsrf_sav(ig,is_oce) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &
2105     &                        - pctsrf_sav(ig,is_sic)
2106        endif
2107      enddo
2108!
2109! Pour rattraper des erreurs d'arrondis
2110!
2111      where (abs(pctsrf_sav(:,is_sic)) .le. 2.*epsilon(pctsrf_sav(1,is_sic)))
2112        pctsrf_sav(:,is_sic) = 0.
2113        pctsrf_sav(:,is_oce) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
2114      endwhere
2115      where (abs(pctsrf_sav(:,is_oce)) .le. 2.*epsilon(pctsrf_sav(1,is_oce)))
2116        pctsrf_sav(:,is_sic) = pctsrf(:,is_oce) + pctsrf(:,is_sic)
2117        pctsrf_sav(:,is_oce) = 0.
2118      endwhere
2119      if (minval(pctsrf_sav(:,is_oce)) < 0.) then
2120        write(*,*)'Pb fraction ocean inferieure a 0'
2121        write(*,*)'au point ',minloc(pctsrf_sav(:,is_oce))
2122        write(*,*)'valeur = ',minval(pctsrf_sav(:,is_oce))
2123        abort_message = 'voir ci-dessus'
2124        call abort_gcm(modname,abort_message,1)
2125      endif
2126      if (minval(pctsrf_sav(:,is_sic)) < 0.) then
2127        write(*,*)'Pb fraction glace inferieure a 0'
2128        write(*,*)'au point ',minloc(pctsrf_sav(:,is_sic))
2129        write(*,*)'valeur = ',minval(pctsrf_sav(:,is_sic))
2130        abort_message = 'voir ci-dessus'
2131        call abort_gcm(modname,abort_message,1)
2132      endif
2133    endif
2134  endif                         ! fin mod(itime, nexca) == 1
2135
2136  if (mod(itime, nexca) == 0) then
2137!
2138! allocation memoire
2139    if (nisurf == is_oce .and. (.not. cumul) ) then
2140      sum_error = 0
2141      allocate(tmp_sols(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2142      allocate(tmp_nsol(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2143      allocate(tmp_rain(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2144      allocate(tmp_snow(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2145      allocate(tmp_evap(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2146      allocate(tmp_tsol(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2147      allocate(tmp_fder(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2148      allocate(tmp_albe(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2149      allocate(tmp_taux(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2150      allocate(tmp_tauy(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2151! -- LOOP
2152       allocate(tmp_windsp(iim,jjphy_nb,2), stat=error); sum_error = sum_error + error
2153! -- LOOP
2154!!$      allocate(tmp_rriv(iim,jjm+1,2), stat=error); sum_error = sum_error + error
2155!!$      allocate(tmp_rcoa(iim,jjm+1,2), stat=error); sum_error = sum_error + error
2156      if (sum_error /= 0) then
2157        abort_message='Pb allocation variables couplees pour l''ecriture'
2158        call abort_gcm(modname,abort_message,1)
2159      endif
2160    endif
2161
2162!
2163! Mise sur la bonne grille des champs a passer au coupleur
2164!
2165    cpl_index = 1
2166    if (nisurf == is_sic) cpl_index = 2
2167    call gath2cpl(cpl_sols(1,cpl_index), tmp_sols(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2168    call gath2cpl(cpl_nsol(1,cpl_index), tmp_nsol(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2169    call gath2cpl(cpl_rain(1,cpl_index), tmp_rain(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2170    call gath2cpl(cpl_snow(1,cpl_index), tmp_snow(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2171    call gath2cpl(cpl_evap(1,cpl_index), tmp_evap(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2172    call gath2cpl(cpl_tsol(1,cpl_index), tmp_tsol(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2173    call gath2cpl(cpl_fder(1,cpl_index), tmp_fder(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2174    call gath2cpl(cpl_albe(1,cpl_index), tmp_albe(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2175    call gath2cpl(cpl_taux(1,cpl_index), tmp_taux(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2176    call gath2cpl(cpl_tauy(1,cpl_index), tmp_tauy(1,1,cpl_index), klon, knon,iim,jjphy_nb,                  knindex)
2177! -- LOOP
2178     call gath2cpl(cpl_windsp(1,cpl_index), tmp_windsp(1,1,cpl_index), klon, knon,iim,jjphy_nb,            knindex)
2179! -- LOOP
2180
2181!
2182! Si le domaine considere est la banquise, on envoie les champs au coupleur
2183!
2184    if (nisurf == is_sic .and. cumul) then
2185      wri_rain = 0.; wri_snow = 0.; wri_rcoa = 0.; wri_rriv = 0.
2186      wri_taux = 0.; wri_tauy = 0.
2187! -- LOOP
2188       wri_windsp = 0.
2189! -- LOOP     
2190      call gath2cpl(pctsrf(1,is_oce), tamp_srf(1,1,1), klon, klon, iim, jjphy_nb, tamp_ind)
2191      call gath2cpl(pctsrf(1,is_sic), tamp_srf(1,1,2), klon, klon, iim, jjphy_nb, tamp_ind)
2192
2193      wri_sol_ice = tmp_sols(:,:,2)
2194      wri_sol_sea = tmp_sols(:,:,1)
2195      wri_nsol_ice = tmp_nsol(:,:,2)
2196      wri_nsol_sea = tmp_nsol(:,:,1)
2197      wri_fder_ice = tmp_fder(:,:,2)
2198      wri_evap_ice = tmp_evap(:,:,2)
2199      wri_evap_sea = tmp_evap(:,:,1)
2200! -- LOOP
2201       wri_windsp = tmp_windsp(:,:,1)
2202! -- LOOP
2203
2204!!$PB
2205      wri_rriv = cpl_rriv(:,:)
2206      wri_rcoa = cpl_rcoa(:,:)
2207
2208!ym  !! ATTENTION ICI
2209     
2210!ym      DO j = 1, jjm + 1
2211!ym        wri_calv(:,j) = sum(cpl_rlic(:,j)) / iim
2212!ym      enddo
2213
2214!Essai OM+JLD : ca marche !!!! (17 mars 2006)
2215      tamp_srf(:,:,3)=0.
2216      CALL gath2cpl( pctsrf(1,is_lic), tamp_srf(1,1,3), klon, klon, iim, jjphy_nb, tamp_ind)
2217
2218!YM pour retrouver resultat avant tamp_srf(:,3)=1.
2219     
2220      DO j = 1, jjphy_nb
2221         wri_calv(:,j) = DOT_PRODUCT (cpl_rlic(1:iim,j), tamp_srf(1:iim,j,3)) / REAL(iim)
2222      ENDDO
2223
2224!ym      wri_calv(:,:)=0.
2225!ym      DO j = 1, jjphy_nb
2226!ym        wri_calv(:,j) = sum(cpl_rlic(:,j))/iim
2227!ym      enddo
2228
2229      IF (.NOT. monocpu) THEN
2230        if (phy_rank /= 0) then
2231#ifdef CPP_PARA
2232          call MPI_RECV(Up,1,MPI_REAL8,phy_rank-1,1234,COMM_LMDZ,status,ierr)
2233          call MPI_SEND(wri_calv(1,1),1,MPI_REAL8,phy_rank-1,1234,COMM_LMDZ,ierr)
2234#endif
2235        endif
2236       
2237        if (phy_rank /= phy_size-1) then
2238#ifdef CPP_PARA
2239          call MPI_SEND(wri_calv(1,jjphy_nb),1,MPI_REAL8,phy_rank+1,1234,COMM_LMDZ,ierr) 
2240          call MPI_RECV(down,1,MPI_REAL8,phy_rank+1,1234,COMM_LMDZ,status,ierr)
2241#endif
2242        endif
2243       
2244        if (phy_rank /=0 .and. iiphy_begin /=1) then
2245          Up=Up+wri_calv(iim,1)
2246          wri_calv(:,1)=Up
2247        endif
2248     
2249        if (phy_rank /=phy_size-1 .and. iiphy_end /= iim) then
2250          Down=Down+wri_calv(1,jjphy_nb)
2251          wri_calv(:,jjphy_nb)=Down     
2252        endif
2253      ENDIF
2254     
2255      where (tamp_zmasq /= 1.)
2256        deno =  tamp_srf(:,:,1) + tamp_srf(:,:,2)
2257        wri_rain = tmp_rain(:,:,1) * tamp_srf(:,:,1) / deno +    &
2258      &            tmp_rain(:,:,2) * tamp_srf(:,:,2) / deno
2259        wri_snow = tmp_snow(:,:,1) * tamp_srf(:,:,1) / deno +    &
2260      &            tmp_snow(:,:,2) * tamp_srf(:,:,2) / deno
2261        wri_taux = tmp_taux(:,:,1) * tamp_srf(:,:,1) / deno +    &
2262      &            tmp_taux(:,:,2) * tamp_srf(:,:,2) / deno
2263        wri_tauy = tmp_tauy(:,:,1) * tamp_srf(:,:,1) / deno +    &
2264      &            tmp_tauy(:,:,2) * tamp_srf(:,:,2) / deno
2265      endwhere
2266!
2267! pour simuler la fonte des glaciers antarctiques
2268!
2269!$$$        wri_rain = wri_rain      &
2270!$$$      &     + coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)
2271!      wri_calv = coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)
2272!
2273! on passe les coordonnees de la grille
2274!
2275
2276!ym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,tmp_lon)
2277!ym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,tmp_lat)
2278
2279      CALL phy2dyn(rlon,tmp_lon,1)
2280      CALL phy2dyn(rlat,tmp_lat,1)
2281     
2282!ym      DO i = 1, iim
2283!ym        tmp_lon(i,1) = rlon(i+1)
2284!ym        tmp_lon(i,jjm + 1) = rlon(i+1)
2285!ym      ENDDO
2286
2287!
2288! sortie netcdf des champs pour le changement de repere
2289!
2290      IF (monocpu) THEN
2291        ndexct(:)=0
2292        CALL histwrite(nidct,'tauxe',itau_w,wri_taux,iim*(jjm+1),ndexct)
2293        CALL histwrite(nidct,'tauyn',itau_w,wri_tauy,iim*(jjm+1),ndexct)
2294        CALL histwrite(nidct,'tmp_lon',itau_w,tmp_lon,iim*(jjm+1),ndexct)
2295        CALL histwrite(nidct,'tmp_lat',itau_w,tmp_lat,iim*(jjm+1),ndexct)
2296      ENDIF 
2297!
2298! calcul 3 coordonnees du vent
2299!
2300      CALL atm2geo (iim , jjphy_nb, wri_taux, wri_tauy, tmp_lon, tmp_lat, &
2301         & wri_tauxx, wri_tauyy, wri_tauzz )
2302!
2303! sortie netcdf des champs apres changement de repere et juste avant
2304! envoi au coupleur
2305!
2306    IF (monocpu) THEN
2307      CALL histwrite(nidct,cl_writ(8),itau_w,wri_sol_ice,iim*(jjm+1),ndexct)
2308      CALL histwrite(nidct,cl_writ(9),itau_w,wri_sol_sea,iim*(jjm+1),ndexct)
2309      CALL histwrite(nidct,cl_writ(10),itau_w,wri_nsol_ice,iim*(jjm+1),ndexct)
2310      CALL histwrite(nidct,cl_writ(11),itau_w,wri_nsol_sea,iim*(jjm+1),ndexct)
2311      CALL histwrite(nidct,cl_writ(12),itau_w,wri_fder_ice,iim*(jjm+1),ndexct)
2312      CALL histwrite(nidct,cl_writ(13),itau_w,wri_evap_ice,iim*(jjm+1),ndexct)
2313      CALL histwrite(nidct,cl_writ(14),itau_w,wri_evap_sea,iim*(jjm+1),ndexct)
2314      CALL histwrite(nidct,cl_writ(15),itau_w,wri_rain,iim*(jjm+1),ndexct)
2315      CALL histwrite(nidct,cl_writ(16),itau_w,wri_snow,iim*(jjm+1),ndexct)
2316      CALL histwrite(nidct,cl_writ(17),itau_w,wri_rcoa,iim*(jjm+1),ndexct)
2317      CALL histwrite(nidct,cl_writ(18),itau_w,wri_rriv,iim*(jjm+1),ndexct)
2318      CALL histwrite(nidct,cl_writ(19),itau_w,wri_calv,iim*(jjm+1),ndexct)
2319      CALL histwrite(nidct,cl_writ(1),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
2320      CALL histwrite(nidct,cl_writ(2),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
2321      CALL histwrite(nidct,cl_writ(3),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
2322      CALL histwrite(nidct,cl_writ(4),itau_w,wri_tauxx,iim*(jjm+1),ndexct)
2323      CALL histwrite(nidct,cl_writ(5),itau_w,wri_tauyy,iim*(jjm+1),ndexct)
2324      CALL histwrite(nidct,cl_writ(6),itau_w,wri_tauzz,iim*(jjm+1),ndexct)
2325! -- LOOP
2326      CALL histwrite(nidct,cl_writ(7),itau_w,wri_windsp,iim*(jjm+1),ndexct)
2327! -- LOOP
2328      CALL histsync(nidct)
2329    ENDIF
2330! pas utile      IF (lafin) CALL histclo(nidct)
2331#ifdef CPP_COUPLE
2332#ifdef CPP_PSMILE
2333      il_time_secs=(itime-1)*dtime
2334     
2335      CALL intocpl(il_time_secs, iim, jjm+1, wri_sol_ice, wri_sol_sea, wri_nsol_ice,&
2336      & wri_nsol_sea, wri_fder_ice, wri_evap_ice, wri_evap_sea, wri_rain, &
2337      & wri_snow, wri_rcoa, wri_rriv, wri_calv, wri_tauxx, wri_tauyy,     &
2338      & wri_tauzz, wri_tauxx, wri_tauyy, wri_tauzz,                       &
2339! -- LOOP
2340      & wri_windsp,lafin)
2341! -- LOOP
2342#else
2343      call intocpl(itime, (jjm+1)*iim, wri_sol_ice, wri_sol_sea, wri_nsol_ice,&
2344      & wri_nsol_sea, wri_fder_ice, wri_evap_ice, wri_evap_sea, wri_rain, &
2345      & wri_snow, wri_rcoa, wri_rriv, wri_calv, wri_tauxx, wri_tauyy,     &
2346      & wri_tauzz, wri_tauxx, wri_tauyy, wri_tauzz,                       &
2347! -- LOOP
2348      & wri_windsp,lafin)
2349! -- LOOP
2350#endif
2351#endif
2352!
2353      cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.
2354      cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.
2355      cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.; cpl_rlic = 0.
2356! -- LOOP
2357      cpl_windsp = 0.
2358! -- LOOP
2359!
2360! deallocation memoire variables temporaires
2361!
2362      sum_error = 0
2363      deallocate(tmp_sols, stat=error); sum_error = sum_error + error
2364      deallocate(tmp_nsol, stat=error); sum_error = sum_error + error
2365      deallocate(tmp_rain, stat=error); sum_error = sum_error + error
2366      deallocate(tmp_snow, stat=error); sum_error = sum_error + error
2367      deallocate(tmp_evap, stat=error); sum_error = sum_error + error
2368      deallocate(tmp_fder, stat=error); sum_error = sum_error + error
2369      deallocate(tmp_tsol, stat=error); sum_error = sum_error + error
2370      deallocate(tmp_albe, stat=error); sum_error = sum_error + error
2371      deallocate(tmp_taux, stat=error); sum_error = sum_error + error
2372      deallocate(tmp_tauy, stat=error); sum_error = sum_error + error
2373! -- LOOP
2374      deallocate(tmp_windsp, stat=error); sum_error = sum_error + error
2375! -- LOOP
2376!!$PB
2377!!$      deallocate(tmp_rriv, stat=error); sum_error = sum_error + error
2378!!$      deallocate(tmp_rcoa, stat=error); sum_error = sum_error + error
2379      if (sum_error /= 0) then
2380        abort_message='Pb deallocation variables couplees'
2381        call abort_gcm(modname,abort_message,1)
2382      endif
2383
2384    endif
2385
2386  endif            ! fin (mod(itime, nexca) == 0)
2387!
2388! on range les variables lues/sauvegardees dans les bonnes variables de sortie
2389!
2390  if (nisurf == is_oce) then
2391    call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjphy_nb, knindex)
2392  else if (nisurf == is_sic) then
2393    call cpl2gath(read_sit, tsurf_new, klon, knon,iim,jjphy_nb, knindex)
2394    call cpl2gath(read_alb_sic, alb_new, klon, knon,iim,jjphy_nb, knindex)
2395  endif
2396  pctsrf_new(:,nisurf) = pctsrf_sav(:,nisurf)
2397 
2398  if (mod(itime, nexca) == -1) then
2399    tmp_field=0.
2400    do i = 1, knon
2401      ig = knindex(i)
2402      tmp_field(ig) = 1.
2403    enddo   
2404!ym    call WriteField_phy('knindex',tmp_field,1)
2405   
2406    tmp_field=0.
2407    do i = 1, knon
2408      ig = knindex(i)
2409      tmp_field(ig) = tsurf_new(i)
2410    enddo   
2411!ym    call WriteField_phy('tsurf_new',tmp_field,1)
2412   
2413    tmp_field=0.
2414    do i = 1, knon
2415      ig = knindex(i)
2416      tmp_field(ig) = alb_new(i)
2417    enddo   
2418!ym    call WriteField_phy('alb_new',tmp_field,1)
2419   
2420!    tmp_field=0.
2421!   do i = 1, knon
2422!      ig = knindex(i)
2423!      tmp_field(ig) = pctsrf_new(i,nisurf)
2424 !   enddo   
2425!ym    call WriteField_phy('pctsrf_new', pctsrf_new(:,nisurf),1)
2426  endif
2427!ym  do j=1,jjphy_nb
2428!ym    do i=1,iim
2429!ym      print *,phy_rank,'read_sst(',i,',',j,')=',read_sst(i,j)
2430!ym    enddo
2431!ym  enddo
2432 
2433!ym  do i=1,knon
2434!ym    print *,phy_rank,'tsurf_new(',i,')=',tsurf_new(i)
2435!ym  enddo
2436!  if (lafin) call quitcpl
2437
2438  END SUBROUTINE interfoce_cpl
2439!
2440!#########################################################################
2441!
2442   SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &
2443 & radsol, fluxo, fluxg, pctsrf, &
2444 & tslab, seaice, pctsrf_slab)
2445!
2446! Cette routine calcule la temperature d'un slab ocean, la glace de mer
2447! et les pourcentages de la maille couverte par l'ocean libre et/ou
2448! la glace de mer pour un "slab" ocean de 50m
2449!
2450! Conception: Laurent Li
2451! Re-ecriture + adaptation LMDZ4: I. Musat
2452!
2453! input:
2454!   klon         nombre total de points de grille
2455!   debut        logical: 1er appel a la physique
2456!   itap         numero du pas de temps
2457!   dtime        pas de temps de la physique (en s)
2458!   ijour        jour dans l'annee en cours
2459!   radsol       rayonnement net au sol (LW + SW)
2460!   fluxo        flux turbulent (sensible + latent) sur les mailles oceaniques
2461!   fluxg        flux de conduction entre la surface de la glace de mer et l'ocean
2462!   pctsrf       tableau des pourcentages de surface de chaque maille
2463! output:
2464!   tslab        temperature de l'ocean libre
2465!   seaice       glace de mer (kg/m2)
2466!   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
2467!
2468#include "indicesol.inc"
2469#include "clesphys.inc"
2470! Parametres d'entree
2471  integer, intent(IN) :: klon
2472  logical, intent(IN) :: debut
2473  INTEGER, intent(IN) :: itap
2474  REAL, intent(IN) :: dtime
2475  INTEGER, intent(IN) :: ijour
2476  REAL, dimension(klon), intent(IN) :: radsol
2477  REAL, dimension(klon), intent(IN) :: fluxo
2478  REAL, dimension(klon), intent(IN) :: fluxg
2479  real, dimension(klon, nbsrf), intent(IN) :: pctsrf
2480! Parametres de sortie
2481  real, dimension(klon), intent(INOUT) :: tslab
2482  real, dimension(klon), intent(INOUT)        :: seaice ! glace de mer (kg/m2)
2483  real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_slab
2484!
2485! Variables locales :
2486  REAL :: amn, amx
2487  INTEGER, save :: lmt_pas, julien, idayvrai
2488!$OMP THREADPRIVATE(lmt_pas, julien, idayvrai)
2489  REAL, parameter :: unjour=86400.
2490  real, allocatable, dimension(:), save :: tmp_tslab, tmp_seaice
2491  REAL, allocatable, dimension(:), save :: slab_bils
2492  REAL, allocatable, dimension(:), save :: lmt_bils
2493!$OMP THREADPRIVATE(tmp_tslab, tmp_seaice,slab_bils,lmt_bils)
2494  logical,save              :: check = .false.
2495!$OMP THREADPRIVATE(check)
2496
2497!
2498  REAL, parameter :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
2499  REAL, parameter :: cbing=0.334e+05        ! J/kg
2500  real, dimension(klon)                 :: siceh !hauteur de la glace de mer (m)
2501  INTEGER :: i
2502  integer :: sum_error, error
2503  REAL :: zz, za, zb
2504!
2505  character (len = 80) :: abort_message
2506  character (len = 20) :: modname = 'interfoce_slab'
2507!
2508  julien = MOD(ijour,360)
2509  sum_error = 0
2510  IF (debut) THEN
2511   allocate(slab_bils(klon), stat = error); sum_error = sum_error + error
2512   allocate(lmt_bils(klon), stat = error); sum_error = sum_error + error
2513   allocate(tmp_tslab(klon), stat = error); sum_error = sum_error + error
2514   allocate(tmp_seaice(klon), stat = error); sum_error = sum_error + error
2515   if (sum_error /= 0) then
2516    abort_message='Pb allocation var. slab_bils,lmt_bils,tmp_tslab,tmp_seaice'
2517    call abort_gcm(modname,abort_message,1)
2518   endif
2519   tmp_tslab=tslab
2520   tmp_seaice=seaice
2521   lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
2522!
2523   IF (check) THEN
2524    PRINT*,'interfoce_slab klon, debut, itap, dtime, ijour, &
2525 &          lmt_pas ', klon, debut, itap, dtime, ijour, &
2526 &          lmt_pas
2527   ENDIF !check
2528!
2529   PRINT*, '************************'
2530   PRINT*, 'SLAB OCEAN est actif, prenez precautions !'
2531   PRINT*, '************************'
2532!
2533! a mettre un slab_bils aussi en force !!!
2534!
2535    DO i = 1, klon
2536     slab_bils(i) = 0.0
2537    ENDDO
2538!
2539   ENDIF !debut
2540!
2541   IF (check ) THEN
2542    amn=MIN(tmp_tslab(1),1000.)
2543    amx=MAX(tmp_tslab(1),-1000.)
2544    DO i=2, klon
2545     amn=MIN(tmp_tslab(i),amn)
2546     amx=MAX(tmp_tslab(i),amx)
2547    ENDDO
2548!
2549    PRINT*,' debut min max tslab',amn,amx
2550!
2551!!
2552    PRINT*,' itap,lmt_pas unjour',itap,lmt_pas,unjour
2553   ENDIF
2554!!
2555!
2556   pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)
2557!
2558! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
2559!
2560   IF (MOD(itap,lmt_pas) .EQ. 1) THEN !1er pas de temps de la journee
2561    idayvrai = ijour
2562    CALL condsurf(julien,idayvrai, lmt_bils)
2563   ENDIF !(MOD(itap-1,lmt_pas) .EQ. 0) THEN
2564
2565   DO i = 1, klon
2566     IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
2567     &  (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
2568!
2569! fabriquer de la glace si congelation atteinte:
2570!
2571      IF (tmp_tslab(i).LT.(RTT-1.8)) THEN
2572       zz =  (RTT-1.8)-tmp_tslab(i)
2573       tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz
2574       seaice(i) = tmp_seaice(i)
2575       tmp_tslab(i) = RTT-1.8
2576      ENDIF
2577!
2578! faire fondre de la glace si temperature est superieure a 0:
2579!
2580      IF ((tmp_tslab(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN
2581       zz = cyang/cbing * (tmp_tslab(i)-RTT)
2582       zz = MIN(zz,tmp_seaice(i))
2583       tmp_seaice(i) = tmp_seaice(i) - zz
2584       seaice(i) = tmp_seaice(i)
2585       tmp_tslab(i) = tmp_tslab(i) - zz*cbing/cyang
2586      ENDIF
2587!
2588! limiter la glace de mer a 10 metres (10000 kg/m2)
2589!
2590      IF(tmp_seaice(i).GT.45.) THEN
2591       tmp_seaice(i) = MIN(tmp_seaice(i),10000.0)
2592      ELSE
2593       tmp_seaice(i) = 0.
2594      ENDIF
2595      seaice(i) = tmp_seaice(i)
2596      siceh(i)=tmp_seaice(i)/1000. !en metres
2597!
2598! determiner la nature du sol (glace de mer ou ocean libre):
2599!
2600! on fait dependre la fraction de seaice "pctsrf(i,is_sic)"
2601! de l'epaisseur de seaice :
2602! pctsrf(i,is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm
2603! et pctsrf(i,is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur
2604!
2605
2606     IF(.NOT.ok_slab_sicOBS) then
2607      pctsrf_slab(i,is_sic)=MIN(siceh(i)/0.20, &
2608     &                      1.-(pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)))
2609      pctsrf_slab(i,is_oce)=1.0 - &
2610     &      (pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)+pctsrf_slab(i,is_sic))
2611     ELSE
2612      IF (i.EQ.1) print*,'cas ok_slab_sicOBS TRUE : passe sur la modif.'
2613     ENDIF !(.NOT.ok_slab_sicOBS) then
2614     ENDIF !pctsrf
2615   ENDDO
2616!
2617! Calculer le bilan du flux de chaleur au sol :
2618!
2619   DO i = 1, klon
2620    za = radsol(i) + fluxo(i)
2621    zb = fluxg(i)
2622    IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
2623   &   (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
2624     slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i,is_oce) &
2625   &             +zb*pctsrf_slab(i,is_sic))/ FLOAT(lmt_pas)
2626    ENDIF
2627   ENDDO !klon
2628!
2629! calcul tslab
2630!
2631   IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee
2632!
2633! calcul tslab
2634    amn=MIN(tmp_tslab(1),1000.)
2635    amx=MAX(tmp_tslab(1),-1000.)
2636    DO i = 1, klon
2637      IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
2638     &    (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
2639       tmp_tslab(i) = tmp_tslab(i) + &
2640     & (slab_bils(i)-lmt_bils(i)) &
2641     &                         /cyang*unjour
2642! on remet l'accumulation a 0
2643       slab_bils(i) = 0.
2644      ENDIF !pctsrf
2645!
2646      IF (check) THEN
2647       IF(i.EQ.1) THEN 
2648        PRINT*,'i tmp_tslab MOD(itap,lmt_pas).EQ.0 sicINTER',i,itap, &
2649      & tmp_tslab(i), tmp_seaice(i)
2650       ENDIF
2651!
2652       amn=MIN(tmp_tslab(i),amn)
2653       amx=MAX(tmp_tslab(i),amx)
2654      ENDIF
2655    ENDDO !klon
2656   ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN
2657!
2658   IF ( check ) THEN
2659    PRINT*,'fin min max tslab',amn,amx
2660   ENDIF
2661!
2662   tslab = tmp_tslab
2663   seaice = tmp_seaice
2664  END SUBROUTINE interfoce_slab
2665!
2666!#########################################################################
2667!
2668  SUBROUTINE interfoce_lim(itime, dtime, jour, &
2669     & klon_xx, nisurf, knon, knindex, &
2670     & debut,  &
2671     & lmt_sst_p, pctsrf_new_p)
2672     
2673     USE dimphy,klon=>klon2,klon2=>klon
2674
2675#include "indicesol.inc"
2676
2677! Cette routine sert d'interface entre le modele atmospherique et un fichier
2678! de conditions aux limites
2679!
2680! L. Fairhead 02/2000
2681!
2682! input:
2683!   itime        numero du pas de temps courant
2684!   dtime        pas de temps de la physique (en s)
2685!   jour         jour a lire dans l'annee
2686!   nisurf       index de la surface a traiter (1 = sol continental)
2687!   knon         nombre de points dans le domaine a traiter
2688!   knindex      index des points de la surface a traiter
2689!   klon         taille de la grille
2690!   debut        logical: 1er appel a la physique (initialisation)
2691!
2692! output:
2693!   lmt_sst      SST lues dans le fichier de CL
2694!   pctsrf_new   sous-maille fractionnelle
2695!
2696
2697
2698! Parametres d'entree
2699  integer, intent(IN) :: itime
2700  real   , intent(IN) :: dtime
2701  integer, intent(IN) :: jour
2702  integer, intent(in) :: klon_xx
2703  integer, intent(IN) :: nisurf
2704  integer, intent(IN) :: knon
2705  integer, dimension(klon2), intent(in) :: knindex
2706  logical, intent(IN) :: debut
2707
2708! Parametres de sortie
2709  real, intent(out), dimension(klon2) :: lmt_sst_p
2710  real, intent(out), dimension(klon2,nbsrf) :: pctsrf_new_p
2711
2712!  real, dimension(klon) :: lmt_sst
2713  real, dimension(klon,nbsrf) :: pctsrf_new
2714
2715! Variables locales
2716  integer     :: ii
2717  INTEGER,save :: lmt_pas     ! frequence de lecture des conditions limites
2718                             ! (en pas de physique)
2719!$OMP THREADPRIVATE(lmt_pas)
2720  logical,save :: deja_lu    ! pour indiquer que le jour a lire a deja
2721                             ! lu pour une surface precedente
2722!$OMP THREADPRIVATE(deja_lu)
2723  integer,save :: jour_lu
2724!$OMP THREADPRIVATE(jour_lu)
2725  integer      :: ierr
2726  character (len = 20) :: modname = 'interfoce_lim'
2727  character (len = 80) :: abort_message
2728  character (len = 20),save :: fich ='limit.nc'
2729!$OMP THREADPRIVATE(fich)
2730  logical, save     :: newlmt = .TRUE.
2731!$OMP THREADPRIVATE(newlmt)
2732  logical, save     :: check = .FALSE.
2733!$OMP THREADPRIVATE(check)
2734! Champs lus dans le fichier de CL
2735  real, allocatable , save, dimension(:) :: sst_lu_p
2736!$OMP THREADPRIVATE(sst_lu_p)
2737  real, allocatable , save, dimension(:) :: sst_lu_mpi
2738
2739  real, allocatable , save, dimension(:,:) :: pct_tmp_p
2740!$OMP THREADPRIVATE(pct_tmp_p)
2741  real, allocatable , save, dimension(:,:) :: pct_tmp_mpi
2742  real, dimension(klon,nbsrf) :: pct_tmp
2743  real, dimension(klon) :: sst_lu
2744  real, dimension(klon) :: nat_lu
2745!
2746! quelques variables pour netcdf
2747!
2748#include "netcdf.inc"
2749  integer              :: nid, nvarid
2750  integer, dimension(2) :: start, epais
2751!
2752! Fin declaration
2753!
2754 
2755  if (debut .and. .not. allocated(sst_lu_p)) then
2756    lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
2757    jour_lu = jour - 1
2758    allocate(sst_lu_p(klon_omp))
2759    allocate(pct_tmp_p(klon_omp,nbsrf))
2760  endif
2761
2762  if ((jour - jour_lu) /= 0) deja_lu = .false.
2763 
2764  if (check) write(*,*)modname,' :: jour, jour_lu, deja_lu', jour, jour_lu, deja_lu
2765  if (check) write(*,*)modname,' :: itime, lmt_pas ', itime, lmt_pas,dtime
2766
2767! Tester d'abord si c'est le moment de lire le fichier
2768  if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then
2769!
2770! Ouverture du fichier
2771!
2772!$OMP MASTER
2773    if (.not. allocated(sst_lu_mpi)) allocate(sst_lu_mpi(klon_mpi))
2774    if (.not. allocated(pct_tmp_mpi)) allocate(pct_tmp_mpi(klon_mpi,nbsrf))
2775   
2776    if (phy_rank==0) then
2777   
2778    fich = trim(fich)
2779    ierr = NF_OPEN (fich, NF_NOWRITE,nid)
2780    if (ierr.NE.NF_NOERR) then
2781      abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
2782      call abort_gcm(modname,abort_message,1)
2783    endif
2784!
2785! La tranche de donnees a lire:
2786!
2787    start(1) = 1
2788    start(2) = jour
2789    epais(1) = klon
2790    epais(2) = 1
2791!
2792    if (newlmt) then
2793!
2794! Fraction "ocean"
2795!
2796      ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)
2797      if (ierr /= NF_NOERR) then
2798        abort_message = 'Le champ <FOCE> est absent'
2799        call abort_gcm(modname,abort_message,1)
2800      endif
2801#ifdef NC_DOUBLE
2802      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_oce))
2803#else
2804      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_oce))
2805#endif
2806      if (ierr /= NF_NOERR) then
2807        abort_message = 'Lecture echouee pour <FOCE>'
2808        call abort_gcm(modname,abort_message,1)
2809      endif
2810!
2811! Fraction "glace de mer"
2812!
2813      ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)
2814      if (ierr /= NF_NOERR) then
2815        abort_message = 'Le champ <FSIC> est absent'
2816        call abort_gcm(modname,abort_message,1)
2817      endif
2818#ifdef NC_DOUBLE
2819      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_sic))
2820#else
2821      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_sic))
2822#endif
2823      if (ierr /= NF_NOERR) then
2824        abort_message = 'Lecture echouee pour <FSIC>'
2825        call abort_gcm(modname,abort_message,1)
2826      endif
2827!
2828! Fraction "terre"
2829!
2830      ierr = NF_INQ_VARID(nid, 'FTER', nvarid)
2831      if (ierr /= NF_NOERR) then
2832        abort_message = 'Le champ <FTER> est absent'
2833        call abort_gcm(modname,abort_message,1)
2834      endif
2835#ifdef NC_DOUBLE
2836      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_ter))
2837#else
2838      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_ter))
2839#endif
2840      if (ierr /= NF_NOERR) then
2841        abort_message = 'Lecture echouee pour <FTER>'
2842        call abort_gcm(modname,abort_message,1)
2843      endif
2844!
2845! Fraction "glacier terre"
2846!
2847      ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)
2848      if (ierr /= NF_NOERR) then
2849        abort_message = 'Le champ <FLIC> est absent'
2850        call abort_gcm(modname,abort_message,1)
2851      endif
2852#ifdef NC_DOUBLE
2853      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais,pct_tmp(1,is_lic))
2854#else
2855      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_lic))
2856#endif
2857      if (ierr /= NF_NOERR) then
2858        abort_message = 'Lecture echouee pour <FLIC>'
2859        call abort_gcm(modname,abort_message,1)
2860      endif
2861!
2862    else  ! on en est toujours a rnatur
2863!
2864      ierr = NF_INQ_VARID(nid, 'NAT', nvarid)
2865      if (ierr /= NF_NOERR) then
2866        abort_message = 'Le champ <NAT> est absent'
2867        call abort_gcm(modname,abort_message,1)
2868      endif
2869#ifdef NC_DOUBLE
2870      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, nat_lu)
2871#else
2872      ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, nat_lu)
2873#endif
2874      if (ierr /= NF_NOERR) then
2875        abort_message = 'Lecture echouee pour <NAT>'
2876        call abort_gcm(modname,abort_message,1)
2877      endif
2878!
2879! Remplissage des fractions de surface
2880! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice
2881!
2882      pct_tmp = 0.0
2883      do ii = 1, klon
2884        pct_tmp(ii,nint(nat_lu(ii)) + 1) = 1.
2885      enddo
2886
2887!
2888!  On se retrouve avec ocean en 1 et terre en 2 alors qu'on veut le contraire
2889!
2890      pctsrf_new = pct_tmp
2891      pctsrf_new (:,2)= pct_tmp (:,1)
2892      pctsrf_new (:,1)= pct_tmp (:,2)
2893      pct_tmp = pctsrf_new
2894    endif ! fin test sur newlmt
2895!
2896! Lecture SST
2897!
2898    ierr = NF_INQ_VARID(nid, 'SST', nvarid)
2899    if (ierr /= NF_NOERR) then
2900      abort_message = 'Le champ <SST> est absent'
2901      call abort_gcm(modname,abort_message,1)
2902    endif
2903#ifdef NC_DOUBLE
2904    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, sst_lu)
2905#else
2906    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, sst_lu)
2907#endif
2908    if (ierr /= NF_NOERR) then
2909      abort_message = 'Lecture echouee pour <SST>'
2910      call abort_gcm(modname,abort_message,1)
2911    endif   
2912
2913!
2914! Fin de lecture
2915!
2916    ierr = NF_CLOSE(nid)
2917   endif ! phyrank
2918!
2919! Recopie des variables dans les champs de sortie
2920!
2921  call ScatterField(sst_lu,sst_lu_mpi,1)
2922  call ScatterField(pct_tmp(:,is_oce),pct_tmp_mpi(:,is_oce),1)
2923  call ScatterField(pct_tmp(:,is_sic),pct_tmp_mpi(:,is_sic),1)
2924!$OMP END MASTER
2925!$OMP BARRIER
2926  call ScatterField_omp(sst_lu_mpi,sst_lu_p,1)
2927  call ScatterField_omp(pct_tmp_mpi(:,is_oce),pct_tmp_p(:,is_oce),1)
2928  call ScatterField_omp(pct_tmp_mpi(:,is_sic),pct_tmp_p(:,is_sic),1)
2929   deja_lu = .true.
2930   jour_lu = jour
2931  endif   
2932 
2933  lmt_sst_p = 999999999.
2934 
2935  do ii = 1, knon
2936    lmt_sst_p(ii) = sst_lu_p(knindex(ii))
2937  enddo
2938
2939  do ii=1,klon2
2940    pctsrf_new_p(ii,is_oce)=pct_tmp_p(ii,is_oce)
2941    pctsrf_new_p(ii,is_sic)=pct_tmp_p(ii,is_sic)
2942  enddo
2943 
2944
2945  END SUBROUTINE interfoce_lim
2946
2947!
2948!#########################################################################
2949!
2950  SUBROUTINE interfsur_lim(itime, dtime, jour, &
2951     & klon_xx, nisurf, knon, knindex, &
2952     & debut,  &
2953     & lmt_alb_p, lmt_rug_p)
2954
2955     USE dimphy,klon=>klon2,klon2=>klon
2956
2957! Cette routine sert d'interface entre le modele atmospherique et un fichier
2958! de conditions aux limites
2959!
2960! L. Fairhead 02/2000
2961!
2962! input:
2963!   itime        numero du pas de temps courant
2964!   dtime        pas de temps de la physique (en s)
2965!   jour         jour a lire dans l'annee
2966!   nisurf       index de la surface a traiter (1 = sol continental)
2967!   knon         nombre de points dans le domaine a traiter
2968!   knindex      index des points de la surface a traiter
2969!   klon         taille de la grille
2970!   debut        logical: 1er appel a la physique (initialisation)
2971!
2972! output:
2973!   lmt_sst      SST lues dans le fichier de CL
2974!   lmt_alb      Albedo lu
2975!   lmt_rug      longueur de rugosite lue
2976!   pctsrf_new   sous-maille fractionnelle
2977!
2978
2979
2980! Parametres d'entree
2981  integer, intent(IN) :: itime
2982  real   , intent(IN) :: dtime
2983  integer, intent(IN) :: jour
2984  integer, intent(IN) :: nisurf
2985  integer, intent(IN) :: knon
2986  integer, intent(IN) :: klon_xx
2987  integer, dimension(klon2), intent(in) :: knindex
2988  logical, intent(IN) :: debut
2989
2990! Parametres de sortie
2991  real, intent(out), dimension(klon2) :: lmt_alb_p
2992  real, intent(out), dimension(klon2) :: lmt_rug_p
2993
2994!  real,  dimension(klon) :: lmt_alb
2995!  real,  dimension(klon) :: lmt_rug
2996
2997! Variables locales
2998  integer     :: ii
2999  integer,save :: lmt_pas     ! frequence de lecture des conditions limites
3000                             ! (en pas de physique)
3001!$OMP THREADPRIVATE(lmt_pas)
3002  logical,save :: deja_lu_sur! pour indiquer que le jour a lire a deja
3003                             ! lu pour une surface precedente
3004!$OMP THREADPRIVATE(deja_lu_sur)
3005  integer,save :: jour_lu_sur
3006!$OMP THREADPRIVATE(jour_lu_sur)
3007  integer      :: ierr
3008  character (len = 20) :: modname = 'interfsur_lim'
3009  character (len = 80) :: abort_message
3010  character (len = 20),save :: fich ='limit.nc'
3011!$OMP THREADPRIVATE(fich)
3012  logical,save     :: newlmt = .false.
3013!$OMP THREADPRIVATE(newlmt)
3014  logical,save     :: check = .false.
3015!$OMP THREADPRIVATE(check)
3016! Champs lus dans le fichier de CL
3017  real, allocatable , save, dimension(:) :: alb_lu_p, rug_lu_p
3018!$OMP THREADPRIVATE(alb_lu_p, rug_lu_p)
3019  real, allocatable , save, dimension(:) :: alb_lu_mpi, rug_lu_mpi
3020  real, dimension(klon) :: alb_lu, rug_lu
3021!
3022! quelques variables pour netcdf
3023!
3024#include "netcdf.inc"
3025  integer ,save             :: nid, nvarid
3026!$OMP THREADPRIVATE(nid, nvarid)
3027  integer, dimension(2),save :: start, epais
3028!$OMP THREADPRIVATE(start, epais)
3029!
3030! Fin declaration
3031!
3032 
3033  if (debut) then
3034    lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
3035    jour_lu_sur = jour - 1
3036    allocate(alb_lu_p(klon_omp))
3037    allocate(rug_lu_p(klon_omp))
3038  endif
3039
3040  if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false.
3041 
3042  if (check) write(*,*)modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, deja_lu_sur
3043  if (check) write(*,*)modname,':: itime, lmt_pas', itime, lmt_pas
3044  if (check) call flush(6)
3045
3046! Tester d'abord si c'est le moment de lire le fichier
3047  if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu_sur) then
3048
3049!$OMP MASTER
3050    if (.not. allocated(alb_lu_mpi)) allocate(alb_lu_mpi(klon_mpi))
3051    if (.not. allocated(rug_lu_mpi)) allocate(rug_lu_mpi(klon_mpi)) 
3052  if (phy_rank==0) then
3053!
3054! Ouverture du fichier
3055!
3056    fich = trim(fich)
3057    IF (check) WRITE(*,*)modname,' ouverture fichier ',fich
3058    if (check) CALL flush(6)
3059    ierr = NF_OPEN (fich, NF_NOWRITE,nid)
3060    if (ierr.NE.NF_NOERR) then
3061      abort_message = 'Pb d''ouverture du fichier de conditions aux limites'
3062      call abort_gcm(modname,abort_message,1)
3063    endif
3064!
3065! La tranche de donnees a lire:
3066 
3067    start(1) = 1
3068    start(2) = jour
3069    epais(1) = klon
3070    epais(2) = 1
3071!
3072! Lecture Albedo
3073!
3074    ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
3075    if (ierr /= NF_NOERR) then
3076      abort_message = 'Le champ <ALB> est absent'
3077      call abort_gcm(modname,abort_message,1)
3078    endif
3079#ifdef NC_DOUBLE
3080    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, alb_lu)
3081#else
3082    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)
3083#endif
3084    if (ierr /= NF_NOERR) then
3085      abort_message = 'Lecture echouee pour <ALB>'
3086      call abort_gcm(modname,abort_message,1)
3087    endif
3088!
3089! Lecture rugosite !
3090    ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
3091    if (ierr /= NF_NOERR) then
3092      abort_message = 'Le champ <RUG> est absent'
3093      call abort_gcm(modname,abort_message,1)
3094    endif
3095#ifdef NC_DOUBLE
3096    ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start,epais, rug_lu)
3097#else
3098    ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)
3099#endif
3100    if (ierr /= NF_NOERR) then
3101      abort_message = 'Lecture echouee pour <RUG>'
3102      call abort_gcm(modname,abort_message,1)
3103    endif
3104
3105!
3106! Fin de lecture
3107!
3108    ierr = NF_CLOSE(nid)
3109
3110
3111  endif  !! phyrank
3112
3113    call ScatterField(alb_lu,alb_lu_mpi,1)
3114    call ScatterField(rug_lu,rug_lu_mpi,1)
3115!$OMP END MASTER
3116!$OMP BARRIER
3117
3118    call ScatterField_omp(alb_lu_mpi,alb_lu_p,1)
3119    call ScatterField_omp(rug_lu_mpi,rug_lu_p,1)
3120   
3121    deja_lu_sur = .true.
3122    jour_lu_sur = jour
3123
3124
3125  endif
3126 
3127!
3128! Recopie des variables dans les champs de sortie
3129!
3130!!$  lmt_alb(:) = 0.0
3131!!$  lmt_rug(:) = 0.0
3132 
3133  lmt_alb_p(:) = 999999.
3134  lmt_rug_p(:) = 999999.
3135  DO ii = 1, knon
3136    lmt_alb_p(ii) = alb_lu_p(knindex(ii))
3137    lmt_rug_p(ii) = rug_lu_p(knindex(ii))
3138  enddo
3139
3140
3141  END SUBROUTINE interfsur_lim
3142
3143!
3144!#########################################################################
3145!
3146
3147  SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime, &
3148     & tsurf, p1lay, cal, beta, coef1lay, ps, &
3149     & precip_rain, precip_snow, snow, qsurf, &
3150     & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
3151     & petAcoef, peqAcoef, petBcoef, peqBcoef, &
3152     & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
3153     USE dimphy,only : omp_rank
3154! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
3155! une temperature de surface (au cas ou ok_veget = false)
3156!
3157! L. Fairhead 4/2000
3158!
3159! input:
3160!   knon         nombre de points a traiter
3161!   nisurf       surface a traiter
3162!   tsurf        temperature de surface
3163!   p1lay        pression 1er niveau (milieu de couche)
3164!   cal          capacite calorifique du sol
3165!   beta         evap reelle
3166!   coef1lay     coefficient d'echange
3167!   ps           pression au sol
3168!   precip_rain  precipitations liquides
3169!   precip_snow  precipitations solides
3170!   snow         champs hauteur de neige
3171!   runoff       runoff en cas de trop plein
3172!   petAcoef     coeff. A de la resolution de la CL pour t
3173!   peqAcoef     coeff. A de la resolution de la CL pour q
3174!   petBcoef     coeff. B de la resolution de la CL pour t
3175!   peqBcoef     coeff. B de la resolution de la CL pour q
3176!   radsol       rayonnement net aus sol (LW + SW)
3177!   dif_grnd     coeff. diffusion vers le sol profond
3178!
3179! output:
3180!   tsurf_new    temperature au sol
3181!   qsurf        humidite de l'air au dessus du sol
3182!   fluxsens     flux de chaleur sensible
3183!   fluxlat      flux de chaleur latente
3184!   dflux_s      derivee du flux de chaleur sensible / Ts
3185!   dflux_l      derivee du flux de chaleur latente  / Ts
3186!
3187
3188#include "YOETHF.inc"
3189#include "FCTTRE.inc"
3190#include "indicesol.inc"
3191#include "YOMCST.inc"
3192
3193! Parametres d'entree
3194  integer, intent(IN) :: knon, nisurf, klon
3195  real   , intent(IN) :: dtime
3196  real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
3197  real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
3198  real, dimension(klon), intent(IN) :: ps, q1lay
3199  real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
3200  real, dimension(klon), intent(IN) :: precip_rain, precip_snow
3201  real, dimension(klon), intent(IN) :: radsol, dif_grnd
3202  real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
3203  real, dimension(klon), intent(INOUT) :: snow, qsurf
3204
3205! Parametres sorties
3206  real, dimension(klon), intent(OUT):: tsurf_new, evap, fluxsens, fluxlat
3207  real, dimension(klon), intent(OUT):: dflux_s, dflux_l
3208
3209! Variables locales
3210  integer :: i
3211  real, dimension(klon) :: zx_mh, zx_nh, zx_oh
3212  real, dimension(klon) :: zx_mq, zx_nq, zx_oq
3213  real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
3214  real, dimension(klon) :: zx_sl, zx_k1
3215  real, dimension(klon) :: zx_q_0 , d_ts
3216  real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
3217  real                  :: bilan_f, fq_fonte
3218  REAL                  :: subli, fsno
3219  REAL                  :: qsat_new, q1_new
3220  real, parameter :: t_grnd = 271.35, t_coup = 273.15
3221!! PB temporaire en attendant mieux pour le modele de neige
3222  REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
3223!
3224  logical, save         :: check = .false.
3225!$OMP THREADPRIVATE(check)
3226  character (len = 20)  :: modname = 'calcul_fluxs'
3227  logical, save         :: fonte_neige = .false.
3228!$OMP THREADPRIVATE(fonte_neige)
3229  real, save            :: max_eau_sol = 150.0
3230!$OMP THREADPRIVATE(max_eau_sol)
3231  character (len = 80) :: abort_message
3232  logical,save         :: first = .true.,second=.false.
3233!$OMP THREADPRIVATE(first,second)
3234
3235  if (check) write(*,*)'Entree ', modname,' surface = ',nisurf
3236
3237  IF (check) THEN
3238      WRITE(*,*)' radsol (min, max)' &
3239         &     , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
3240      CALL flush(6)
3241  ENDIF
3242
3243  if (size(coastalflow) /= knon .AND. nisurf == is_ter) then
3244    write(*,*)'Bizarre, le nombre de points continentaux'
3245    write(*,*)'a change entre deux appels. J''arrete ...'
3246    abort_message='Pb run_off'
3247    call abort_gcm(modname,abort_message,1)
3248  endif
3249!
3250! Traitement neige et humidite du sol
3251!
3252!!$  WRITE(*,*)'test calcul_flux, surface ', nisurf
3253!!PB test
3254!!$    if (nisurf == is_oce) then
3255!!$      snow = 0.
3256!!$      qsol = max_eau_sol
3257!!$    else
3258!!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)
3259!!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))
3260!!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)
3261!!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime
3262!!$    endif
3263!!$    IF (nisurf /= is_ter) qsol = max_eau_sol
3264
3265
3266!
3267! Initialisation
3268!
3269  evap = 0.
3270  fluxsens=0.
3271  fluxlat=0.
3272  dflux_s = 0.
3273  dflux_l = 0. 
3274!
3275! zx_qs = qsat en kg/kg
3276!
3277  DO i = 1, knon
3278    zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
3279  IF (thermcep) THEN
3280      zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
3281      zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
3282      zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
3283      zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
3284      zx_qs=MIN(0.5,zx_qs)
3285      zcor=1./(1.-retv*zx_qs)
3286      zx_qs=zx_qs*zcor
3287      zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
3288     &                 /RLVTT / zx_pkh(i)
3289    ELSE
3290      IF (tsurf(i).LT.t_coup) THEN
3291        zx_qs = qsats(tsurf(i)) / ps(i)
3292        zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
3293     &                    / zx_pkh(i)
3294      ELSE
3295        zx_qs = qsatl(tsurf(i)) / ps(i)
3296        zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
3297     &               / zx_pkh(i)
3298      ENDIF
3299    ENDIF
3300    zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
3301    zx_qsat(i) = zx_qs
3302    zx_coef(i) = coef1lay(i) &
3303     & * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
3304     & * p1lay(i)/(RD*t1lay(i))
3305
3306  ENDDO
3307
3308
3309! === Calcul de la temperature de surface ===
3310!
3311! zx_sl = chaleur latente d'evaporation ou de sublimation
3312!
3313  do i = 1, knon
3314    zx_sl(i) = RLVTT
3315    if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
3316    zx_k1(i) = zx_coef(i)
3317  enddo
3318
3319
3320  do i = 1, knon
3321! Q
3322    zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
3323    zx_mq(i) = beta(i) * zx_k1(i) * &
3324     &             (peqAcoef(i) - zx_qsat(i) &
3325     &                          + zx_dq_s_dt(i) * tsurf(i)) &
3326     &             / zx_oq(i)
3327    zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
3328     &                              / zx_oq(i)
3329
3330! H
3331    zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
3332    zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
3333    zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
3334
3335! Tsurface
3336    tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
3337     &             (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
3338     &                 + dif_grnd(i) * t_grnd * dtime)/ &
3339     &          ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
3340     &                       zx_nh(i) + zx_sl(i) * zx_nq(i)) & 
3341     &                     + dtime * dif_grnd(i))
3342
3343!
3344! Y'a-t-il fonte de neige?
3345!
3346!    fonte_neige = (nisurf /= is_oce) .AND. &
3347!     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
3348!     & .AND. (tsurf_new(i) >= RTT)
3349!    if (fonte_neige) tsurf_new(i) = RTT 
3350    d_ts(i) = tsurf_new(i) - tsurf(i)
3351!    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
3352!    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
3353!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
3354!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
3355    evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
3356    fluxlat(i) = - evap(i) * zx_sl(i)
3357    fluxsens(i) = zx_mh(i) + zx_nh(i) * tsurf_new(i)
3358! Derives des flux dF/dTs (W m-2 K-1):
3359    dflux_s(i) = zx_nh(i)
3360    dflux_l(i) = (zx_sl(i) * zx_nq(i))
3361! Nouvelle valeure de l'humidite au dessus du sol
3362    qsat_new=zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
3363    q1_new = peqAcoef(i) - peqBcoef(i)*evap(i)*dtime
3364    qsurf(i)=q1_new*(1.-beta(i)) + beta(i)*qsat_new
3365!
3366! en cas de fonte de neige
3367!
3368!    if (fonte_neige) then
3369!      bilan_f = radsol(i) + fluxsens(i) - (zx_sl(i) * evap (i)) - &
3370!     &          dif_grnd(i) * (tsurf_new(i) - t_grnd) - &
3371!     &          RCPD * (zx_pkh(i))/cal(i)/dtime * (tsurf_new(i) - tsurf(i))
3372!      bilan_f = max(0., bilan_f)
3373!      fq_fonte = bilan_f / zx_sl(i)
3374!      snow(i) = max(0., snow(i) - fq_fonte * dtime)
3375!      qsol(i) = qsol(i) + (fq_fonte * dtime)
3376!    endif
3377!!$    if (nisurf == is_ter)  &
3378!!$     &  run_off(i) = run_off(i) + max(qsol(i) - max_eau_sol, 0.0)
3379!!$    qsol(i) = min(qsol(i), max_eau_sol)
3380  ENDDO
3381
3382  END SUBROUTINE calcul_fluxs
3383!
3384!#########################################################################
3385!
3386  SUBROUTINE gath2cpl(champ_in, champ_out, klon, knon, iim, jmp1, knindex)
3387  use dimphy, only: liste_i,liste_j,jjphy_begin,jjphy_nb,phy_rank,phy_size
3388  implicit none
3389 
3390! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
3391! au coupleur.
3392!
3393!
3394! input:         
3395!   champ_in     champ sur la grille gathere       
3396!   knon         nombre de points dans le domaine a traiter
3397!   knindex      index des points de la surface a traiter
3398!   klon         taille de la grille
3399!   iim,jjm      dimension de la grille 2D
3400!
3401! output:
3402!   champ_out    champ sur la grille 2D
3403!
3404! input
3405  integer                   :: klon, knon, iim, jmp1
3406  real, dimension(klon)     :: champ_in
3407  integer, dimension(klon)  :: knindex
3408! output
3409  real, dimension(iim,jmp1)  :: champ_out
3410! local
3411  integer                   :: i, ig, j
3412  real, dimension(klon)     :: tamp
3413
3414  tamp = 0.
3415  do i = 1, knon
3416    ig = knindex(i)
3417    tamp(ig) = champ_in(i)
3418  enddo   
3419
3420!ym  ig = 1
3421!ym  champ_out(:,1) = tamp(ig)
3422!ym  do j = 2, jjm
3423!ym    do i = 1, iim
3424!ym      ig = ig + 1
3425!ym      champ_out(i,j) = tamp(ig)
3426!ym    enddo
3427!ym  enddo
3428!ym  ig = ig + 1
3429!ym  champ_out(:,jjm+1) = tamp(ig)
3430
3431  do ig=1,klon
3432    i=liste_i(ig)
3433    j=liste_j(ig)-jjphy_begin+1
3434    champ_out(i,j)=tamp(ig)
3435  enddo
3436 
3437  if (phy_rank==0) champ_out(:,1)=tamp(1)
3438  if (phy_rank==phy_size-1) champ_out(:,jjphy_nb)=tamp(klon)
3439
3440  END SUBROUTINE gath2cpl
3441!
3442!#########################################################################
3443!
3444  SUBROUTINE cpl2gath(champ_in, champ_out, klon, knon, iim, jmp1, knindex)
3445  use dimphy, only : liste_i, liste_j, jjphy_begin
3446  implicit none
3447! Cette routine ecrit un champ 'gathered' sur la grille 2D pour le passer
3448! au coupleur.
3449!
3450!
3451! input:         
3452!   champ_in     champ sur la grille gathere       
3453!   knon         nombre de points dans le domaine a traiter
3454!   knindex      index des points de la surface a traiter
3455!   klon         taille de la grille
3456!   iim,jjm      dimension de la grille 2D
3457!
3458! output:
3459!   champ_out    champ sur la grille 2D
3460!
3461! input
3462  integer                   :: klon, knon, iim, jmp1
3463  real, dimension(iim,jmp1)     :: champ_in
3464  integer, dimension(klon)  :: knindex
3465! output
3466  real, dimension(klon)  :: champ_out
3467! local
3468  integer                   :: i, ig, j
3469  real, dimension(klon)     :: tamp
3470  logical ,save                  :: check = .false.
3471!$OMP THREADPRIVATE(check)
3472
3473!ym  ig = 1
3474!ym  tamp(ig) = champ_in(1,1)
3475!ym  do j = 2, jjm
3476!ym    do i = 1, iim
3477!ym      ig = ig + 1
3478!ym      tamp(ig) = champ_in(i,j)
3479!ym    enddo
3480!ym  enddo
3481!ym  ig = ig + 1
3482!ym  tamp(ig) = champ_in(1,jjm+1)
3483
3484  do ig=1,klon
3485   i=liste_i(ig)
3486   j=liste_j(ig)-jjphy_begin+1
3487   tamp(ig)=champ_in(i,j)
3488  enddo
3489 
3490  do i = 1, knon
3491    ig = knindex(i)
3492    champ_out(i) = tamp(ig)
3493  enddo   
3494
3495  END SUBROUTINE cpl2gath
3496!
3497!#########################################################################
3498!
3499  SUBROUTINE albsno(klon, knon,dtime,agesno,alb_neig_grid, precip_snow)
3500  IMPLICIT none
3501 
3502  INTEGER :: klon, knon
3503  INTEGER, PARAMETER :: nvm = 8
3504  REAL   :: dtime
3505  REAL, dimension(klon,nvm) :: veget
3506  REAL, DIMENSION(klon) :: alb_neig_grid, agesno, precip_snow
3507 
3508  INTEGER :: i, nv
3509 
3510  REAL, DIMENSION(nvm),SAVE :: init, decay
3511!$OMP THREADPRIVATE(init, decay)
3512  REAL :: as
3513  DATA init /0.55, 0.14, 0.18, 0.29, 0.15, 0.15, 0.14, 0./
3514  DATA decay/0.30, 0.67, 0.63, 0.45, 0.40, 0.14, 0.06, 1./
3515 
3516  veget = 0.
3517  veget(:,1) = 1.     ! desert partout
3518  DO i = 1, knon
3519    alb_neig_grid(i) = 0.0
3520  ENDDO
3521  DO nv = 1, nvm
3522    DO i = 1, knon
3523      as = init(nv)+decay(nv)*EXP(-agesno(i)/5.)
3524      alb_neig_grid(i) = alb_neig_grid(i) + veget(i,nv)*as
3525    ENDDO
3526  ENDDO
3527!
3528!! modilation en fonction de l'age de la neige
3529!
3530  DO i = 1, knon
3531    agesno(i)  = (agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
3532            &             * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/0.3)
3533    agesno(i) =  MAX(agesno(i),0.0)
3534  ENDDO
3535 
3536  END SUBROUTINE albsno
3537!
3538!#########################################################################
3539!
3540
3541  SUBROUTINE fonte_neige( klon, knon, nisurf, dtime, &
3542     & tsurf, p1lay, cal, beta, coef1lay, ps, &
3543     & precip_rain, precip_snow, snow, qsol, &
3544     & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
3545     & petAcoef, peqAcoef, petBcoef, peqBcoef, &
3546     & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
3547     & fqcalving,fqfonte,ffonte,run_off_lic_0)
3548
3549! Routine de traitement de la fonte de la neige dans le cas du traitement
3550! de sol simplifie !
3551! LF 03/2001
3552! input:
3553!   knon         nombre de points a traiter
3554!   nisurf       surface a traiter
3555!   tsurf        temperature de surface
3556!   p1lay        pression 1er niveau (milieu de couche)
3557!   cal          capacite calorifique du sol
3558!   beta         evap reelle
3559!   coef1lay     coefficient d'echange
3560!   ps           pression au sol
3561!   precip_rain  precipitations liquides
3562!   precip_snow  precipitations solides
3563!   snow         champs hauteur de neige
3564!   qsol         hauteur d'eau contenu dans le sol
3565!   runoff       runoff en cas de trop plein
3566!   petAcoef     coeff. A de la resolution de la CL pour t
3567!   peqAcoef     coeff. A de la resolution de la CL pour q
3568!   petBcoef     coeff. B de la resolution de la CL pour t
3569!   peqBcoef     coeff. B de la resolution de la CL pour q
3570!   radsol       rayonnement net aus sol (LW + SW)
3571!   dif_grnd     coeff. diffusion vers le sol profond
3572!
3573! output:
3574!   tsurf_new    temperature au sol
3575!   fluxsens     flux de chaleur sensible
3576!   fluxlat      flux de chaleur latente
3577!   dflux_s      derivee du flux de chaleur sensible / Ts
3578!   dflux_l      derivee du flux de chaleur latente  / Ts
3579! in/out:
3580!   run_off_lic_0 run off glacier du pas de temps precedent
3581!
3582
3583#include "YOETHF.inc"
3584!rv#include "FCTTRE.inc"
3585#include "indicesol.inc"
3586!IM cf JLD
3587#include "YOMCST.inc"
3588
3589! Parametres d'entree
3590  integer, intent(IN) :: knon, nisurf, klon
3591  real   , intent(IN) :: dtime
3592  real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
3593  real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
3594  real, dimension(klon), intent(IN) :: ps, q1lay
3595  real, dimension(klon), intent(IN) :: tsurf, p1lay, cal, beta, coef1lay
3596  real, dimension(klon), intent(IN) :: precip_rain, precip_snow
3597  real, dimension(klon), intent(IN) :: radsol, dif_grnd
3598  real, dimension(klon), intent(IN) :: t1lay, u1lay, v1lay
3599  real, dimension(klon), intent(INOUT) :: snow, qsol
3600
3601! Parametres sorties
3602  real, dimension(klon), intent(INOUT):: tsurf_new, evap, fluxsens, fluxlat
3603  real, dimension(klon), intent(INOUT):: dflux_s, dflux_l
3604! Flux thermique utiliser pour fondre la neige
3605  real, dimension(klon), intent(INOUT):: ffonte
3606! Flux d'eau "perdu" par la surface et necessaire pour que limiter la
3607! hauteur de neige, en kg/m2/s. Et flux d'eau de fonte de la calotte.
3608  REAL, DIMENSION(klon), INTENT(INOUT):: fqcalving, fqfonte
3609  real, dimension(klon), intent(INOUT):: run_off_lic_0
3610! Variables locales
3611! Masse maximum de neige (kg/m2). Au dessus de ce seuil, la neige
3612! en exces "s'ecoule" (calving)
3613!  real, parameter :: snow_max=1.
3614!IM cf JLD/GK
3615  real, parameter :: snow_max=3000.
3616  integer :: i
3617  real, dimension(klon) :: zx_mh, zx_nh, zx_oh
3618  real, dimension(klon) :: zx_mq, zx_nq, zx_oq
3619  real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
3620  real, dimension(klon) :: zx_sl, zx_k1
3621  real, dimension(klon) :: zx_q_0 , d_ts
3622  real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
3623  real                  :: bilan_f, fq_fonte
3624  REAL                  :: subli, fsno
3625  REAL, DIMENSION(klon) :: bil_eau_s, snow_evap
3626  real, parameter :: t_grnd = 271.35, t_coup = 273.15
3627!! PB temporaire en attendant mieux pour le modele de neige
3628! REAL, parameter :: chasno = RLMLT/(2.3867E+06*0.15)
3629  REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
3630!IM cf JLD/ GKtest
3631  REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)
3632! fin GKtest
3633!
3634  logical, save         :: check = .FALSE.
3635!$OMP THREADPRIVATE(check)
3636  character (len = 20)  :: modname = 'fonte_neige'
3637  logical, save         :: neige_fond = .false.
3638!$OMP THREADPRIVATE(neige_fond)
3639  real, save            :: max_eau_sol = 150.0
3640!$OMP THREADPRIVATE(max_eau_sol)
3641  character (len = 80) :: abort_message
3642  logical,save         :: first = .true.,second=.false.
3643!$OMP THREADPRIVATE(first,second)
3644  real                 :: coeff_rel
3645#include "FCTTRE.inc"
3646
3647
3648  if (check) write(*,*)'Entree ', modname,' surface = ',nisurf
3649
3650! Initialisations
3651  coeff_rel = dtime/(tau_calv * rday)
3652  bil_eau_s(:) = 0.
3653  DO i = 1, knon
3654    zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
3655  IF (thermcep) THEN
3656      zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))
3657      zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
3658      zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
3659      zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)
3660      zx_qs=MIN(0.5,zx_qs)
3661      zcor=1./(1.-retv*zx_qs)
3662      zx_qs=zx_qs*zcor
3663      zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &
3664     &                 /RLVTT / zx_pkh(i)
3665    ELSE
3666      IF (tsurf(i).LT.t_coup) THEN
3667        zx_qs = qsats(tsurf(i)) / ps(i)
3668        zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &
3669     &                    / zx_pkh(i)
3670      ELSE
3671        zx_qs = qsatl(tsurf(i)) / ps(i)
3672        zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &
3673     &               / zx_pkh(i)
3674      ENDIF
3675    ENDIF
3676    zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
3677    zx_qsat(i) = zx_qs
3678    zx_coef(i) = coef1lay(i) &
3679     & * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
3680     & * p1lay(i)/(RD*t1lay(i))
3681  ENDDO
3682
3683
3684! === Calcul de la temperature de surface ===
3685!
3686! zx_sl = chaleur latente d'evaporation ou de sublimation
3687!
3688  do i = 1, knon
3689    zx_sl(i) = RLVTT
3690    if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
3691    zx_k1(i) = zx_coef(i)
3692  enddo
3693
3694
3695  do i = 1, knon
3696! Q
3697    zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
3698    zx_mq(i) = beta(i) * zx_k1(i) * &
3699     &             (peqAcoef(i) - zx_qsat(i) &
3700     &                          + zx_dq_s_dt(i) * tsurf(i)) &
3701     &             / zx_oq(i)
3702    zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
3703     &                              / zx_oq(i)
3704
3705! H
3706    zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
3707    zx_mh(i) = zx_k1(i) * petAcoef(i) / zx_oh(i)
3708    zx_nh(i) = - (zx_k1(i) * RCPD * zx_pkh(i))/ zx_oh(i)
3709  enddo
3710
3711
3712  WHERE (precip_snow > 0.) snow = snow + (precip_snow * dtime)
3713  snow_evap = 0.
3714  WHERE (evap > 0. )
3715    snow_evap = MIN (snow / dtime, evap)
3716    snow = snow - snow_evap * dtime
3717    snow = MAX(0.0, snow)
3718  end where
3719 
3720!  bil_eau_s = bil_eau_s + (precip_rain * dtime) - (evap - snow_evap) * dtime
3721  bil_eau_s = (precip_rain * dtime) - (evap - snow_evap) * dtime
3722
3723!
3724! Y'a-t-il fonte de neige?
3725!
3726  ffonte=0.
3727  do i = 1, knon
3728    neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
3729     & .AND. tsurf_new(i) >= RTT)
3730    if (neige_fond) then
3731      fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i))
3732      ffonte(i) = fq_fonte * RLMLT/dtime
3733      fqfonte(i) = fq_fonte/dtime
3734      snow(i) = max(0., snow(i) - fq_fonte)
3735      bil_eau_s(i) = bil_eau_s(i) + fq_fonte
3736      tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno 
3737!IM cf JLD OK     
3738!IM cf JLD/ GKtest fonte aussi pour la glace
3739      IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
3740        fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,0.0)
3741        ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime
3742        IF ( ok_lic_melt ) THEN
3743           fqfonte(i) = fqfonte(i) + fq_fonte/dtime
3744           bil_eau_s(i) = bil_eau_s(i) + fq_fonte
3745        ENDIF
3746        tsurf_new(i) = RTT
3747      ENDIF
3748      d_ts(i) = tsurf_new(i) - tsurf(i)
3749    endif
3750!
3751!   s'il y a une hauteur trop importante de neige, elle s'coule
3752    fqcalving(i) = max(0., snow(i) - snow_max)/dtime
3753    snow(i)=min(snow(i),snow_max)
3754!
3755    IF (nisurf == is_ter) then
3756      qsol(i) = qsol(i) + bil_eau_s(i)
3757      run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)
3758      qsol(i) = MIN(qsol(i), max_eau_sol)
3759    else if (nisurf == is_lic) then
3760      run_off_lic(i) = (coeff_rel *  fqcalving(i)) + &
3761 &                        (1. - coeff_rel) * run_off_lic_0(i)
3762      run_off_lic_0(i) = run_off_lic(i)
3763      run_off_lic(i) = run_off_lic(i) + fqfonte(i)/dtime
3764    endif
3765  enddo
3766
3767  END SUBROUTINE fonte_neige
3768!
3769!#########################################################################
3770!
3771  END MODULE interface_surf
Note: See TracBrowser for help on using the repository browser.