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

Last change on this file since 223 was 223, checked in by lmdzadmin, 23 years ago

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