/[lmdze]/trunk/phylmd/Interface_surf/interface_surf.f
ViewVC logotype

Diff of /trunk/phylmd/Interface_surf/interface_surf.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC revision 53 by guez, Fri Oct 7 13:11:58 2011 UTC
# Line 1  Line 1 
1  MODULE interface_surf  MODULE interface_surf
2    
3    ! From phylmd/interface_surf.F90,v 1.8 2005/05/25 13:10:09    ! From phylmd/interface_surf.F90, version 1.8 2005/05/25 13:10:09
4    
5    ! Ce module regroupe toutes les routines gérant l'interface entre le modèle    ! Ce module regroupe toutes les routines gérant l'interface entre le
6    ! atmosphérique et les modèles de surface (sols continentaux,    ! modèle atmosphérique et les modèles de surface (sols continentaux,
7    ! océans, glaces).    ! océans, glaces). Les routines sont les suivantes. "interfsurf_hq" :
8    ! Les routines sont les suivantes:    ! routine d'aiguillage vers les interfaces avec les différents
9    ! interfsurf_hq : routine d'aiguillage vers les interfaces avec les    ! modèles de surface ; "interfoce_*" : routines d'interface proprement
10    !    différents modèles de surface    ! dites.
   ! interfoce_* : routines d'interface proprement dites  
11    
12    ! L. Fairhead, LMD, 02/2000    ! L. Fairhead, LMD, february 2000
13    
14    IMPLICIT none    IMPLICIT none
15    
16    PRIVATE    PRIVATE
17    PUBLIC :: interfsurf_hq    PUBLIC :: interfsurf_hq
18    
19    ! run_off      ruissellement total    ! run_off ruissellement total
20    REAL, ALLOCATABLE, DIMENSION(:),SAVE    :: run_off, run_off_lic    REAL, ALLOCATABLE, DIMENSION(:), SAVE :: run_off, run_off_lic
21    real, allocatable, dimension(:),save    :: coastalflow, riverflow    real, allocatable, dimension(:), save :: coastalflow, riverflow
22    
23    REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: tmp_rriv, tmp_rcoa,tmp_rlic    REAL, ALLOCATABLE, DIMENSION(:, :), SAVE :: tmp_rriv, tmp_rcoa, tmp_rlic
24    !! pour simuler la fonte des glaciers antarctiques    ! pour simuler la fonte des glaciers antarctiques
25    REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: coeff_iceberg    REAL, save :: tau_calv
   real, save                              :: surf_maille  
   real, save                              :: cte_flux_iceberg = 6.3e7  
   integer, save                           :: num_antarctic = 1  
   REAL, save                              :: tau_calv  
26    
27  CONTAINS  CONTAINS
28    
29    SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, &    SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, klon, iim, jjm, &
30         & klon, iim, jjm, nisurf, knon, knindex, pctsrf, &         nisurf, knon, knindex, pctsrf, rlon, rlat, cufi, cvfi, debut, lafin, &
31         & rlon, rlat, cufi, cvfi,&         ok_veget, soil_model, nsoilmx, tsoil, qsol, zlev, u1_lay, v1_lay, &
32         & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, qsol,&         temp_air, spechum, epot_air, ccanopy, tq_cdrag, petAcoef, peqAcoef, &
33         & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &         petBcoef, peqBcoef, precip_rain, precip_snow, sollw, sollwdown, swnet, &
34         & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &         swdown, fder, taux, tauy, windsp, rugos, rugoro, albedo, snow, qsurf, &
35         & precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &         tsurf, p1lay, ps, radsol, ocean, npas, nexca, zmasq, evap, fluxsens, &
36         & fder, taux, tauy, &         fluxlat, dflux_l, dflux_s, tsol_rad, tsurf_new, alb_new, alblw, &
37         & windsp, &         emis_new, z0_new, pctsrf_new, agesno, fqcalving, ffonte, &
38         & rugos, rugoro, &         run_off_lic_0, flux_o, flux_g, tslab, seaice)
        & albedo, snow, qsurf, &  
        & tsurf, p1lay, ps, radsol, &  
        & ocean, npas, nexca, zmasq, &  
        & evap, fluxsens, fluxlat, dflux_l, dflux_s, &                
        & tsol_rad, tsurf_new, alb_new, alblw, emis_new, &  
        & z0_new, pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0,&  
        !IM "slab" ocean  
        & flux_o, flux_g, tslab, seaice)  
39    
40      ! Cette routine sert d'aiguillage entre l'atmosphère et la surface      ! Cette routine sert d'aiguillage entre l'atmosphère et la surface
41      ! en général (sols continentaux, océans, glaces) pour les fluxs de      ! en général (sols continentaux, océans, glaces) pour les flux de
42      ! chaleur et d'humidité.      ! chaleur et d'humidité.
43      ! En pratique l'interface se fait entre la couche limite du modèle      ! En pratique l'interface se fait entre la couche limite du modèle
44      ! atmosphérique ("clmain.F") et les routines de surface      ! atmosphérique ("clmain.F") et les routines de surface
# Line 59  CONTAINS Line 46  CONTAINS
46    
47      ! L.Fairhead 02/2000      ! L.Fairhead 02/2000
48    
     ! input:  
     !   klon         nombre total de points de grille  
     !   iim, jjm     nbres de pts de grille  
     !   dtime        pas de temps de la physique (en s)  
     !   date0        jour initial  
     !   jour         jour dans l'annee en cours,  
     !   rmu0         cosinus de l'angle solaire zenithal  
     !   nexca        pas de temps couplage  
     !   nisurf       index de la surface a traiter (1 = sol continental)  
     !   knon         nombre de points de la surface a traiter  
     !   knindex      index des points de la surface a traiter  
     !   pctsrf       tableau des pourcentages de surface de chaque maille  
     !   rlon         longitudes  
     !   rlat         latitudes  
     !   cufi,cvfi    resolution des mailles en x et y (m)  
     !   debut        logical: 1er appel a la physique  
     !   lafin        logical: dernier appel a la physique  
     !   ok_veget     logical: appel ou non au schema de surface continental  
     !                  (si false calcul simplifie des fluxs sur les continents)  
     !   zlev         hauteur de la premiere couche  
     !   u1_lay       vitesse u 1ere couche  
     !   v1_lay       vitesse v 1ere couche  
     !   temp_air     temperature de l'air 1ere couche  
     !   spechum      humidite specifique 1ere couche  
     !   epot_air     temp potentielle de l'air  
     !   ccanopy      concentration CO2 canopee  
     !   tq_cdrag     cdrag  
     !   petAcoef     coeff. A de la resolution de la CL pour t  
     !   peqAcoef     coeff. A de la resolution de la CL pour q  
     !   petBcoef     coeff. B de la resolution de la CL pour t  
     !   peqBcoef     coeff. B de la resolution de la CL pour q  
     !   precip_rain  precipitation liquide  
     !   precip_snow  precipitation solide  
     !   sollw        flux IR net a la surface  
     !   sollwdown    flux IR descendant a la surface  
     !   swnet        flux solaire net  
     !   swdown       flux solaire entrant a la surface  
     !   albedo       albedo de la surface  
     !   tsurf        temperature de surface  
     !   tslab        temperature slab ocean  
     !   pctsrf_slab  pourcentages (0-1) des sous-surfaces dans le slab  
     !   tmp_pctsrf_slab = pctsrf_slab  
     !   p1lay        pression 1er niveau (milieu de couche)  
     !   ps           pression au sol  
     !   radsol       rayonnement net aus sol (LW + SW)  
     !   ocean        type d'ocean utilise (force, slab, couple)  
     !   fder         derivee des flux (pour le couplage)  
     !   taux, tauy   tension de vents  
     !   windsp       module du vent a 10m  
     !   rugos        rugosite  
     !   zmasq        masque terre/ocean  
     !   rugoro       rugosite orographique  
     !   run_off_lic_0 runoff glacier du pas de temps precedent  
   
     ! output:  
     !   evap         evaporation totale  
     !   fluxsens     flux de chaleur sensible  
     !   fluxlat      flux de chaleur latente  
     !   tsol_rad      
     !   tsurf_new    temperature au sol  
     !   alb_new      albedo  
     !   emis_new     emissivite  
     !   z0_new       surface roughness  
     !   pctsrf_new   nouvelle repartition des surfaces  
   
49      use abort_gcm_m, only: abort_gcm      use abort_gcm_m, only: abort_gcm
50      use gath_cpl, only: gath2cpl      use gath_cpl, only: gath2cpl
51      use indicesol      use indicesol
52      use YOMCST      use SUPHEC_M
53      use albsno_m, only: albsno      use albsno_m, only: albsno
54    
55      ! Parametres d'entree      ! Parametres d'entree
56      integer, intent(IN) :: itime     !  numero du pas de temps      ! input:
57        ! klon nombre total de points de grille
58        ! iim, jjm nbres de pts de grille
59        ! dtime pas de temps de la physique (en s)
60        ! date0 jour initial
61        ! jour jour dans l'annee en cours,
62        ! rmu0 cosinus de l'angle solaire zenithal
63        ! nexca pas de temps couplage
64        ! nisurf index de la surface a traiter (1 = sol continental)
65        ! knon nombre de points de la surface a traiter
66        ! knindex index des points de la surface a traiter
67        ! pctsrf tableau des pourcentages de surface de chaque maille
68        ! rlon longitudes
69        ! rlat latitudes
70        ! cufi, cvfi resolution des mailles en x et y (m)
71        ! debut logical: 1er appel a la physique
72        ! lafin logical: dernier appel a la physique
73        ! ok_veget logical: appel ou non au schema de surface continental
74        ! (si false calcul simplifie des fluxs sur les continents)
75        ! zlev hauteur de la premiere couche
76        ! u1_lay vitesse u 1ere couche
77        ! v1_lay vitesse v 1ere couche
78        ! temp_air temperature de l'air 1ere couche
79        ! spechum humidite specifique 1ere couche
80        ! epot_air temp potentielle de l'air
81        ! ccanopy concentration CO2 canopee
82        ! tq_cdrag cdrag
83        ! petAcoef coeff. A de la resolution de la CL pour t
84        ! peqAcoef coeff. A de la resolution de la CL pour q
85        ! petBcoef coeff. B de la resolution de la CL pour t
86        ! peqBcoef coeff. B de la resolution de la CL pour q
87        ! precip_rain precipitation liquide
88        ! precip_snow precipitation solide
89        ! sollw flux IR net a la surface
90        ! sollwdown flux IR descendant a la surface
91        ! swnet flux solaire net
92        ! swdown flux solaire entrant a la surface
93        ! albedo albedo de la surface
94        ! tsurf temperature de surface
95        ! tslab temperature slab ocean
96        ! pctsrf_slab pourcentages (0-1) des sous-surfaces dans le slab
97        ! tmp_pctsrf_slab = pctsrf_slab
98        ! p1lay pression 1er niveau (milieu de couche)
99        ! ps pression au sol
100        ! radsol rayonnement net aus sol (LW + SW)
101        ! ocean type d'ocean utilise ("force" ou "slab" mais pas "couple")
102        ! fder derivee des flux (pour le couplage)
103        ! taux, tauy tension de vents
104        ! windsp module du vent a 10m
105        ! rugos rugosite
106        ! zmasq masque terre/ocean
107        ! rugoro rugosite orographique
108        ! run_off_lic_0 runoff glacier du pas de temps precedent
109        integer, intent(IN) :: itime ! numero du pas de temps
110      integer, intent(IN) :: iim, jjm      integer, intent(IN) :: iim, jjm
111      integer, intent(IN) :: klon      integer, intent(IN) :: klon
112      real, intent(IN) :: dtime      real, intent(IN) :: dtime
113      real, intent(IN) :: date0      real, intent(IN) :: date0
114      integer, intent(IN) :: jour      integer, intent(IN) :: jour
115      real, intent(IN)    :: rmu0(klon)      real, intent(IN) :: rmu0(klon)
116      integer, intent(IN) :: nisurf      integer, intent(IN) :: nisurf
117      integer, intent(IN) :: knon      integer, intent(IN) :: knon
118      integer, dimension(klon), intent(in) :: knindex      integer, dimension(klon), intent(in) :: knindex
119      real, dimension(klon,nbsrf), intent(IN) :: pctsrf      real, dimension(klon, nbsrf), intent(IN) :: pctsrf
120      logical, intent(IN) :: debut, lafin, ok_veget      logical, intent(IN) :: debut, lafin, ok_veget
121      real, dimension(klon), intent(IN) :: rlon, rlat      real, dimension(klon), intent(IN) :: rlon, rlat
122      real, dimension(klon), intent(IN) :: cufi, cvfi      real, dimension(klon), intent(IN) :: cufi, cvfi
# Line 160  CONTAINS Line 135  CONTAINS
135      real, dimension(klon), intent(INOUT) :: tslab      real, dimension(klon), intent(INOUT) :: tslab
136      real, allocatable, dimension(:), save :: tmp_tslab      real, allocatable, dimension(:), save :: tmp_tslab
137      real, dimension(klon), intent(OUT) :: flux_o, flux_g      real, dimension(klon), intent(OUT) :: flux_o, flux_g
138      real, dimension(klon), intent(INOUT)      :: seaice ! glace de mer (kg/m2)      real, dimension(klon), intent(INOUT) :: seaice ! glace de mer (kg/m2)
139      REAL, DIMENSION(klon), INTENT(INOUT) :: radsol,fder      REAL, DIMENSION(klon), INTENT(INOUT) :: radsol, fder
140      real, dimension(klon), intent(IN) :: zmasq      real, dimension(klon), intent(IN) :: zmasq
141      real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro      real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro
142      real, dimension(klon), intent(IN) :: windsp      real, dimension(klon), intent(IN) :: windsp
143      character (len = 6)  :: ocean      character(len=*), intent(IN):: ocean
144      integer              :: npas, nexca ! nombre et pas de temps couplage      integer :: npas, nexca ! nombre et pas de temps couplage
145      real, dimension(klon), intent(INOUT) :: evap, snow, qsurf      real, dimension(klon), intent(INOUT) :: evap, snow, qsurf
146      !! PB ajout pour soil      !! PB ajout pour soil
147      logical          :: soil_model      logical, intent(in):: soil_model
148      integer          :: nsoilmx      integer :: nsoilmx
149      REAL, DIMENSION(klon, nsoilmx) :: tsoil      REAL, DIMENSION(klon, nsoilmx) :: tsoil
150      REAL, dimension(klon), intent(INOUT) :: qsol      REAL, dimension(klon), intent(INOUT) :: qsol
151      REAL, dimension(klon)          :: soilcap      REAL, dimension(klon) :: soilcap
152      REAL, dimension(klon)          :: soilflux      REAL, dimension(klon) :: soilflux
153    
154      ! Parametres de sortie      ! Parametres de sortie
155        ! output:
156        ! evap evaporation totale
157        ! fluxsens flux de chaleur sensible
158        ! fluxlat flux de chaleur latente
159        ! tsol_rad
160        ! tsurf_new temperature au sol
161        ! alb_new albedo
162        ! emis_new emissivite
163        ! z0_new surface roughness
164        ! pctsrf_new nouvelle repartition des surfaces
165      real, dimension(klon), intent(OUT):: fluxsens, fluxlat      real, dimension(klon), intent(OUT):: fluxsens, fluxlat
166      real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new      real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new
167      real, dimension(klon), intent(OUT):: alblw      real, dimension(klon), intent(OUT):: alblw
168      real, dimension(klon), intent(OUT):: emis_new, z0_new      real, dimension(klon), intent(OUT):: emis_new, z0_new
169      real, dimension(klon), intent(OUT):: dflux_l, dflux_s      real, dimension(klon), intent(OUT):: dflux_l, dflux_s
170      real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new      real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_new
171      real, dimension(klon), intent(INOUT):: agesno      real, dimension(klon), intent(INOUT):: agesno
172      real, dimension(klon), intent(INOUT):: run_off_lic_0      real, dimension(klon), intent(INOUT):: run_off_lic_0
173    
174      ! Flux thermique utiliser pour fondre la neige      ! Flux thermique utiliser pour fondre la neige
175      !jld a rajouter   real, dimension(klon), intent(INOUT):: ffonte      !jld a rajouter real, dimension(klon), intent(INOUT):: ffonte
176      real, dimension(klon), intent(INOUT):: ffonte      real, dimension(klon), intent(INOUT):: ffonte
177      ! Flux d'eau "perdue" par la surface et nécessaire pour que limiter la      ! Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
178      ! hauteur de neige, en kg/m2/s      ! hauteur de neige, en kg/m2/s
179      !jld a rajouter   real, dimension(klon), intent(INOUT):: fqcalving      !jld a rajouter real, dimension(klon), intent(INOUT):: fqcalving
180      real, dimension(klon), intent(INOUT):: fqcalving      real, dimension(klon), intent(INOUT):: fqcalving
181      !IM: "slab" ocean - Local      !IM: "slab" ocean - Local
182      real, parameter :: t_grnd=271.35      real, parameter :: t_grnd=271.35
# Line 198  CONTAINS Line 184  CONTAINS
184      integer i      integer i
185      real, allocatable, dimension(:), save :: tmp_flux_o, tmp_flux_g      real, allocatable, dimension(:), save :: tmp_flux_o, tmp_flux_g
186      real, allocatable, dimension(:), save :: tmp_radsol      real, allocatable, dimension(:), save :: tmp_radsol
187      real, allocatable, dimension(:,:), save :: tmp_pctsrf_slab      real, allocatable, dimension(:, :), save :: tmp_pctsrf_slab
188      real, allocatable, dimension(:), save :: tmp_seaice      real, allocatable, dimension(:), save :: tmp_seaice
189    
190      ! Local      ! Local
191      character (len = 20),save :: modname = 'interfsurf_hq'      character (len = 20), save :: modname = 'interfsurf_hq'
192      character (len = 80) :: abort_message      character (len = 80) :: abort_message
193      logical, save        :: first_call = .true.      logical, save :: first_call = .true.
194      integer, save        :: error      integer, save :: error
195      integer              :: ii      integer :: ii
196      logical,save              :: check = .false.      logical, save :: check = .false.
197      real, dimension(klon):: cal, beta, dif_grnd, capsol      real, dimension(klon):: cal, beta, dif_grnd, capsol
198  !!$PB  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.      real, parameter :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.
199      real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.      real, parameter :: calsno=1./(2.3867e+06*.15)
     real, parameter      :: calsno=1./(2.3867e+06*.15)  
200      real, dimension(klon):: tsurf_temp      real, dimension(klon):: tsurf_temp
201      real, dimension(klon):: alb_neig, alb_eau      real, dimension(klon):: alb_neig, alb_eau
202      real, DIMENSION(klon):: zfra      real, DIMENSION(klon):: zfra
203      logical              :: cumul = .false.      logical :: cumul = .false.
204      INTEGER,dimension(1) :: iloc      INTEGER, dimension(1) :: iloc
205      real, dimension(klon):: fder_prev      real, dimension(klon):: fder_prev
206      REAL, dimension(klon) :: bidule      REAL, dimension(klon) :: bidule
207    
208      !-------------------------------------------------------------      !-------------------------------------------------------------
209    
210      if (check) write(*,*) 'Entree ', modname      if (check) write(*, *) 'Entree ', modname
211    
212      ! On doit commencer par appeler les schemas de surfaces continentales      ! On doit commencer par appeler les schemas de surfaces continentales
213      ! car l'ocean a besoin du ruissellement qui est y calcule      ! car l'ocean a besoin du ruissellement qui est y calcule
# Line 230  CONTAINS Line 215  CONTAINS
215      if (first_call) then      if (first_call) then
216         call conf_interface(tau_calv)         call conf_interface(tau_calv)
217         if (nisurf /= is_ter .and. klon > 1) then         if (nisurf /= is_ter .and. klon > 1) then
218            write(*,*)' *** Warning ***'            write(*, *)' *** Warning ***'
219            write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter            write(*, *)' nisurf = ', nisurf, ' /= is_ter = ', is_ter
220            write(*,*)'or on doit commencer par les surfaces continentales'            write(*, *)'or on doit commencer par les surfaces continentales'
221            abort_message='voir ci-dessus'            abort_message='voir ci-dessus'
222            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
223         endif         endif
224         if (ocean /= 'slab' .and. ocean /= 'force' .and. ocean /= 'couple') then         if (ocean /= 'slab' .and. ocean /= 'force') then
225            write(*,*)' *** Warning ***'            write(*, *)' *** Warning ***'
226            write(*,*)'Option couplage pour l''ocean = ', ocean            write(*, *)'Option couplage pour l''ocean = ', ocean
227            abort_message='option pour l''ocean non valable'            abort_message='option pour l''ocean non valable'
228            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
229         endif         endif
230         if ( is_oce > is_sic ) then         if ( is_oce > is_sic ) then
231            write(*,*)' *** Warning ***'            write(*, *)' *** Warning ***'
232            write(*,*)' Pour des raisons de sequencement dans le code'            write(*, *)' Pour des raisons de sequencement dans le code'
233            write(*,*)' l''ocean doit etre traite avant la banquise'            write(*, *)' l''ocean doit etre traite avant la banquise'
234            write(*,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic            write(*, *)' or is_oce = ', is_oce, '> is_sic = ', is_sic
235            abort_message='voir ci-dessus'            abort_message='voir ci-dessus'
236            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
237         endif         endif
238      endif      endif
239      first_call = .false.      first_call = .false.
240    
241      ! Initialisations diverses      ! Initialisations diverses
242      !  
243      ffonte(1:knon)=0.      ffonte(1:knon)=0.
244      fqcalving(1:knon)=0.      fqcalving(1:knon)=0.
245    
# Line 266  CONTAINS Line 251  CONTAINS
251      !IM: "slab" ocean; initialisations      !IM: "slab" ocean; initialisations
252      flux_o = 0.      flux_o = 0.
253      flux_g = 0.      flux_g = 0.
254      !  
255      if (.not. allocated(tmp_flux_o)) then      if (.not. allocated(tmp_flux_o)) then
256         allocate(tmp_flux_o(klon), stat = error)         allocate(tmp_flux_o(klon), stat = error)
257         DO i=1, knon         DO i=1, knon
# Line 274  CONTAINS Line 259  CONTAINS
259         ENDDO         ENDDO
260         if (error /= 0) then         if (error /= 0) then
261            abort_message='Pb allocation tmp_flux_o'            abort_message='Pb allocation tmp_flux_o'
262            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
263         endif         endif
264      endif      endif
265      if (.not. allocated(tmp_flux_g)) then        if (.not. allocated(tmp_flux_g)) then
266         allocate(tmp_flux_g(klon), stat = error)         allocate(tmp_flux_g(klon), stat = error)
267         DO i=1, knon         DO i=1, knon
268            tmp_flux_g(knindex(i))=flux_g(i)            tmp_flux_g(knindex(i))=flux_g(i)
269         ENDDO         ENDDO
270         if (error /= 0) then         if (error /= 0) then
271            abort_message='Pb allocation tmp_flux_g'            abort_message='Pb allocation tmp_flux_g'
272            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
273         endif         endif
274      endif      endif
275      if (.not. allocated(tmp_radsol)) then        if (.not. allocated(tmp_radsol)) then
276         allocate(tmp_radsol(klon), stat = error)         allocate(tmp_radsol(klon), stat = error)
277         if (error /= 0) then         if (error /= 0) then
278            abort_message='Pb allocation tmp_radsol'            abort_message='Pb allocation tmp_radsol'
279            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
280         endif         endif
281      endif      endif
282      DO i=1, knon      DO i=1, knon
283         tmp_radsol(knindex(i))=radsol(i)         tmp_radsol(knindex(i))=radsol(i)
284      ENDDO      ENDDO
285      if (.not. allocated(tmp_pctsrf_slab)) then      if (.not. allocated(tmp_pctsrf_slab)) then
286         allocate(tmp_pctsrf_slab(klon,nbsrf), stat = error)         allocate(tmp_pctsrf_slab(klon, nbsrf), stat = error)
287         if (error /= 0) then         if (error /= 0) then
288            abort_message='Pb allocation tmp_pctsrf_slab'            abort_message='Pb allocation tmp_pctsrf_slab'
289            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
290         endif         endif
291         DO i=1, klon         DO i=1, klon
292            tmp_pctsrf_slab(i,1:nbsrf)=pctsrf(i,1:nbsrf)            tmp_pctsrf_slab(i, 1:nbsrf)=pctsrf(i, 1:nbsrf)
293         ENDDO         ENDDO
294      endif      endif
295      !  
296      if (.not. allocated(tmp_seaice)) then      if (.not. allocated(tmp_seaice)) then
297         allocate(tmp_seaice(klon), stat = error)         allocate(tmp_seaice(klon), stat = error)
298         if (error /= 0) then         if (error /= 0) then
299            abort_message='Pb allocation tmp_seaice'            abort_message='Pb allocation tmp_seaice'
300            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
301         endif         endif
302         DO i=1, klon         DO i=1, klon
303            tmp_seaice(i)=seaice(i)            tmp_seaice(i)=seaice(i)
304         ENDDO         ENDDO
305      endif      endif
306      !  
307      if (.not. allocated(tmp_tslab)) then      if (.not. allocated(tmp_tslab)) then
308         allocate(tmp_tslab(klon), stat = error)         allocate(tmp_tslab(klon), stat = error)
309         if (error /= 0) then         if (error /= 0) then
310            abort_message='Pb allocation tmp_tslab'            abort_message='Pb allocation tmp_tslab'
311            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
312         endif         endif
313      endif      endif
314      DO i=1, klon      DO i=1, klon
315         tmp_tslab(i)=tslab(i)         tmp_tslab(i)=tslab(i)
316      ENDDO      ENDDO
317      !  
318      ! Aiguillage vers les differents schemas de surface      ! Aiguillage vers les differents schemas de surface
319    
320      if (nisurf == is_ter) then      if (nisurf == is_ter) then
321         !  
322         ! Surface "terre" appel a l'interface avec les sols continentaux         ! Surface "terre" appel a l'interface avec les sols continentaux
323         !  
324         ! allocation du run-off         ! allocation du run-off
325         if (.not. allocated(coastalflow)) then         if (.not. allocated(coastalflow)) then
326            allocate(coastalflow(knon), stat = error)            allocate(coastalflow(knon), stat = error)
327            if (error /= 0) then            if (error /= 0) then
328               abort_message='Pb allocation coastalflow'               abort_message='Pb allocation coastalflow'
329               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
330            endif            endif
331            allocate(riverflow(knon), stat = error)            allocate(riverflow(knon), stat = error)
332            if (error /= 0) then            if (error /= 0) then
333               abort_message='Pb allocation riverflow'               abort_message='Pb allocation riverflow'
334               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
335            endif            endif
336            allocate(run_off(knon), stat = error)            allocate(run_off(knon), stat = error)
337            if (error /= 0) then            if (error /= 0) then
338               abort_message='Pb allocation run_off'               abort_message='Pb allocation run_off'
339               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
340            endif            endif
341            !cym                  !cym
342            run_off=0.0            run_off=0.0
343            !cym            !cym
344    
345  !!$PB  !!$PB
346            ALLOCATE (tmp_rriv(iim,jjm+1), stat=error)            ALLOCATE (tmp_rriv(iim, jjm+1), stat=error)
347            if (error /= 0) then            if (error /= 0) then
348               abort_message='Pb allocation tmp_rriv'               abort_message='Pb allocation tmp_rriv'
349               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
350            endif            endif
351            ALLOCATE (tmp_rcoa(iim,jjm+1), stat=error)            ALLOCATE (tmp_rcoa(iim, jjm+1), stat=error)
352            if (error /= 0) then            if (error /= 0) then
353               abort_message='Pb allocation tmp_rcoa'               abort_message='Pb allocation tmp_rcoa'
354               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
355            endif            endif
356            ALLOCATE (tmp_rlic(iim,jjm+1), stat=error)            ALLOCATE (tmp_rlic(iim, jjm+1), stat=error)
357            if (error /= 0) then            if (error /= 0) then
358               abort_message='Pb allocation tmp_rlic'               abort_message='Pb allocation tmp_rlic'
359               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
360            endif            endif
361            tmp_rriv = 0.0            tmp_rriv = 0.0
362            tmp_rcoa = 0.0            tmp_rcoa = 0.0
# Line 379  CONTAINS Line 364  CONTAINS
364    
365  !!$  !!$
366         else if (size(coastalflow) /= knon) then         else if (size(coastalflow) /= knon) then
367            write(*,*)'Bizarre, le nombre de points continentaux'            write(*, *)'Bizarre, le nombre de points continentaux'
368            write(*,*)'a change entre deux appels. J''arrete ...'            write(*, *)'a change entre deux appels. J''arrete ...'
369            abort_message='voir ci-dessus'            abort_message='voir ci-dessus'
370            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
371         endif         endif
372         coastalflow = 0.         coastalflow = 0.
373         riverflow = 0.         riverflow = 0.
374         !  
375         ! Calcul age de la neige         ! Calcul age de la neige
376    
377         if (.not. ok_veget) then         if (.not. ok_veget) then
378            !            ! calcul albedo: lecture albedo fichier boundary conditions
379            ! calcul albedo: lecture albedo fichier CL puis ajout albedo neige            ! puis ajout albedo neige
380            !            call interfsur_lim(itime, dtime, jour, klon, nisurf, knon, knindex, &
381            call interfsur_lim(itime, dtime, jour, &                 debut, alb_new, z0_new)
382                 & klon, nisurf, knon, knindex, debut,  &  
                & alb_new, z0_new)  
           !    
383            ! calcul snow et qsurf, hydrol adapté            ! calcul snow et qsurf, hydrol adapté
           !  
384            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
385    
386            IF (soil_model) THEN            IF (soil_model) THEN
387               CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)               CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, &
388                      soilflux)
389               cal(1:knon) = RCPD / soilcap(1:knon)               cal(1:knon) = RCPD / soilcap(1:knon)
390               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)               radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
391            ELSE            ELSE
392               cal = RCPD * capsol               cal = RCPD * capsol
 !!$      cal = capsol  
393            ENDIF            ENDIF
394            CALL calcul_fluxs( klon, knon, nisurf, dtime, &            CALL calcul_fluxs( klon, knon, nisurf, dtime, &
395                 &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &                 tsurf, p1lay, cal, beta, tq_cdrag, ps, &
396                 &   precip_rain, precip_snow, snow, qsurf,  &                 precip_rain, precip_snow, snow, qsurf, &
397                 &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &                 radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
398                 &   petAcoef, peqAcoef, petBcoef, peqBcoef, &                 petAcoef, peqAcoef, petBcoef, peqBcoef, &
399                 &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)                 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
400    
401            CALL fonte_neige( klon, knon, nisurf, dtime, &            CALL fonte_neige( klon, knon, nisurf, dtime, &
402                 &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &                 tsurf, p1lay, cal, beta, tq_cdrag, ps, &
403                 &   precip_rain, precip_snow, snow, qsol,  &                 precip_rain, precip_snow, snow, qsol, &
404                 &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &                 radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
405                 &   petAcoef, peqAcoef, petBcoef, peqBcoef, &                 petAcoef, peqAcoef, petBcoef, peqBcoef, &
406                 &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &                 tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
407                 &   fqcalving,ffonte, run_off_lic_0)                 fqcalving, ffonte, run_off_lic_0)
408    
409            call albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))              call albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)
410            where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.            where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
411            zfra(1:knon) = max(0.0,min(1.0,snow(1:knon)/(snow(1:knon)+10.0)))            zfra(1:knon) = max(0.0, min(1.0, snow(1:knon)/(snow(1:knon)+10.0)))
412            alb_new(1 : knon)  = alb_neig(1 : knon) *zfra(1:knon) + &            alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
413                 &                     alb_new(1 : knon)*(1.0-zfra(1:knon))                 alb_new(1 : knon)*(1.0-zfra(1:knon))
414            z0_new = sqrt(z0_new**2+rugoro**2)            z0_new = sqrt(z0_new**2+rugoro**2)
415            alblw(1 : knon) = alb_new(1 : knon)            alblw(1 : knon) = alb_new(1 : knon)
   
        else  
416         endif         endif
        !  
        ! Remplissage des pourcentages de surface  
        !  
        pctsrf_new(:,nisurf) = pctsrf(:,nisurf)  
417    
418           ! Remplissage des pourcentages de surface
419           pctsrf_new(:, nisurf) = pctsrf(:, nisurf)
420      else if (nisurf == is_oce) then      else if (nisurf == is_oce) then
   
        if (check) write(*,*)'ocean, nisurf = ',nisurf  
   
        !  
421         ! Surface "ocean" appel a l'interface avec l'ocean         ! Surface "ocean" appel a l'interface avec l'ocean
422         !         if (ocean == 'slab') then
        if (ocean == 'couple') then  
           if (nexca == 0) then  
              abort_message='nexca = 0 dans interfoce_cpl'  
              call abort_gcm(modname,abort_message,1)  
           endif  
   
           cumul = .false.  
   
           iloc = maxloc(fder(1:klon))  
           if (check) then  
              if (fder(iloc(1))> 0.) then  
                 WRITE(*,*)'**** Debug fder ****'  
                 WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))  
              endif  
           endif  
 !!$  
 !!$      where(fder.gt.0.)  
 !!$        fder = 0.  
 !!$      endwhere  
   
           call interfoce_cpl(itime, dtime, cumul, &  
                & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &  
                & ocean, npas, nexca, debut, lafin, &  
                & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &  
                & fluxlat, fluxsens, fder, albedo, taux, tauy, &  
                & windsp, &  
                & zmasq, &  
                & tsurf_new, alb_new, &  
                & pctsrf_new)  
   
           !IM: "slab" ocean  
        else if (ocean == 'slab  ') then  
423            tsurf_new = tsurf            tsurf_new = tsurf
424            pctsrf_new = tmp_pctsrf_slab            pctsrf_new = tmp_pctsrf_slab
425            !         else
426         else                              ! lecture conditions limites            ! lecture conditions limites
427            call interfoce_lim(itime, dtime, jour, &            call interfoce_lim(itime, dtime, jour, klon, nisurf, knon, knindex, &
428                 &  klon, nisurf, knon, knindex, &                 debut, tsurf_new, pctsrf_new)
                &  debut, &  
                &  tsurf_new, pctsrf_new)  
   
429         endif         endif
430    
431         tsurf_temp = tsurf_new         tsurf_temp = tsurf_new
432         cal = 0.         cal = 0.
433         beta = 1.         beta = 1.
434         dif_grnd = 0.         dif_grnd = 0.
435         alb_neig(:) = 0.         alb_neig = 0.
436         agesno(:) = 0.         agesno = 0.
437    
438         call calcul_fluxs( klon, knon, nisurf, dtime, &         call calcul_fluxs( klon, knon, nisurf, dtime, &
439              &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &              tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
440              &   precip_rain, precip_snow, snow, qsurf,  &              precip_rain, precip_snow, snow, qsurf, &
441              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
442              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              petAcoef, peqAcoef, petBcoef, peqBcoef, &
443              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)              tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
444    
445         fder_prev = fder             fder_prev = fder
446         fder = fder_prev + dflux_s + dflux_l         fder = fder_prev + dflux_s + dflux_l
447    
448         iloc = maxloc(fder(1:klon))         iloc = maxloc(fder(1:klon))
449         if (check.and.fder(iloc(1))> 0.) then         if (check.and.fder(iloc(1))> 0.) then
450            WRITE(*,*)'**** Debug fder****'            WRITE(*, *)'**** Debug fder****'
451            WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))            WRITE(*, *)'max fder(', iloc(1), ') = ', fder(iloc(1))
452            WRITE(*,*)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &            WRITE(*, *)'fder_prev, dflux_s, dflux_l', fder_prev(iloc(1)), &
453                 &                        dflux_s(iloc(1)), dflux_l(iloc(1))                 dflux_s(iloc(1)), dflux_l(iloc(1))
454         endif         endif
 !!$  
 !!$      where(fder.gt.0.)  
 !!$        fder = 0.  
 !!$      endwhere  
455    
456         !IM: flux ocean-atmosphere utile pour le "slab" ocean         !IM: flux ocean-atmosphere utile pour le "slab" ocean
457         DO i=1, knon         DO i=1, knon
# Line 526  CONTAINS Line 461  CONTAINS
461            tmp_flux_o(knindex(i)) = flux_o(i)            tmp_flux_o(knindex(i)) = flux_o(i)
462            tmp_radsol(knindex(i))=radsol(i)            tmp_radsol(knindex(i))=radsol(i)
463         ENDDO         ENDDO
464         !  
465         ! 2eme appel a interfoce pour le cumul des champs (en particulier         ! 2eme appel a interfoce pour le cumul des champs (en particulier
466         ! fluxsens et fluxlat calcules dans calcul_fluxs)         ! fluxsens et fluxlat calcules dans calcul_fluxs)
        !  
        if (ocean == 'couple') then  
   
           cumul = .true.  
467    
468            call interfoce_cpl(itime, dtime, cumul, &         if (ocean == 'slab ') then
                & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &  
                & ocean, npas, nexca, debut, lafin, &  
                & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &  
                & fluxlat, fluxsens, fder, albedo, taux, tauy, &  
                & windsp, &  
                & zmasq, &  
                & tsurf_new, alb_new, &  
                & pctsrf_new)  
   
           !IM: "slab" ocean  
        else if (ocean == 'slab  ') then  
           !  
469            seaice=tmp_seaice            seaice=tmp_seaice
470            cumul = .true.            cumul = .true.
471            call interfoce_slab(klon, debut, itime, dtime, jour, &            call interfoce_slab(klon, debut, itime, dtime, jour, &
472                 & tmp_radsol, tmp_flux_o, tmp_flux_g, pctsrf, &                 tmp_radsol, tmp_flux_o, tmp_flux_g, pctsrf, &
473                 & tslab, seaice, pctsrf_new)                 tslab, seaice, pctsrf_new)
474            !  
475            tmp_pctsrf_slab=pctsrf_new            tmp_pctsrf_slab=pctsrf_new
476            DO i=1, knon            DO i=1, knon
477               tsurf_new(i)=tslab(knindex(i))               tsurf_new(i)=tslab(knindex(i))
478            ENDDO !i            ENDDO
           !  
479         endif         endif
480    
        !  
481         ! calcul albedo         ! calcul albedo
        !  
   
482         if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then         if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then
483            CALL alboc(FLOAT(jour),rlat,alb_eau)            CALL alboc(FLOAT(jour), rlat, alb_eau)
484         else  ! cycle diurne         else ! cycle diurne
485            CALL alboc_cd(rmu0,alb_eau)            CALL alboc_cd(rmu0, alb_eau)
486         endif         endif
487         DO ii =1, knon         DO ii =1, knon
488            alb_new(ii) = alb_eau(knindex(ii))            alb_new(ii) = alb_eau(knindex(ii))
# Line 575  CONTAINS Line 490  CONTAINS
490    
491         z0_new = sqrt(rugos**2 + rugoro**2)         z0_new = sqrt(rugos**2 + rugoro**2)
492         alblw(1:knon) = alb_new(1:knon)         alblw(1:knon) = alb_new(1:knon)
   
        !  
493      else if (nisurf == is_sic) then      else if (nisurf == is_sic) then
494           if (check) write(*, *)'sea ice, nisurf = ', nisurf
495    
        if (check) write(*,*)'sea ice, nisurf = ',nisurf  
   
        !  
496         ! Surface "glace de mer" appel a l'interface avec l'ocean         ! Surface "glace de mer" appel a l'interface avec l'ocean
        !  
        !  
        if (ocean == 'couple') then  
   
           cumul =.false.  
   
           iloc = maxloc(fder(1:klon))  
           if (check.and.fder(iloc(1))> 0.) then  
              WRITE(*,*)'**** Debug fder ****'  
              WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))  
           endif  
 !!$  
 !!$      where(fder.gt.0.)  
 !!$        fder = 0.  
 !!$      endwhere  
   
           call interfoce_cpl(itime, dtime, cumul, &  
                & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &  
                & ocean, npas, nexca, debut, lafin, &  
                & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &  
                & fluxlat, fluxsens, fder, albedo, taux, tauy, &  
                & windsp, &  
                & zmasq, &  
                & tsurf_new, alb_new, &  
                & pctsrf_new)  
497    
           tsurf_temp = tsurf_new  
           cal = 0.  
           dif_grnd = 0.  
           beta = 1.0  
498    
499            !IM: "slab" ocean         if (ocean == 'slab ') then
        else if (ocean == 'slab  ') then  
500            pctsrf_new=tmp_pctsrf_slab            pctsrf_new=tmp_pctsrf_slab
501            !  
502            DO ii = 1, knon            DO ii = 1, knon
503               tsurf_new(ii) = tsurf(ii)               tsurf_new(ii) = tsurf(ii)
504               IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then               IF (pctsrf_new(knindex(ii), nisurf) < EPSFRA) then
505                  snow(ii) = 0.0                  snow(ii) = 0.0
506                  tsurf_new(ii) = RTT - 1.8                  tsurf_new(ii) = RTT - 1.8
507                  IF (soil_model) tsoil(ii,:) = RTT -1.8                  IF (soil_model) tsoil(ii, :) = RTT -1.8
508               ENDIF               ENDIF
509            ENDDO            ENDDO
510    
511            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
512    
513            IF (soil_model) THEN            IF (soil_model) THEN
514               CALL soil(dtime, nisurf, knon,snow, tsurf_new, tsoil,soilcap, soilflux)               CALL soil(dtime, nisurf, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
515               cal(1:knon) = RCPD / soilcap(1:knon)               cal(1:knon) = RCPD / soilcap(1:knon)
516               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)               radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
517            ELSE            ELSE
518               dif_grnd = 1.0 / tau_gl               dif_grnd = 1.0 / tau_gl
519               cal = RCPD * calice               cal = RCPD * calice
# Line 640  CONTAINS Line 521  CONTAINS
521            ENDIF            ENDIF
522            tsurf_temp = tsurf_new            tsurf_temp = tsurf_new
523            beta = 1.0            beta = 1.0
524            !  
525         ELSE         ELSE
526            !                              ! lecture conditions limites            ! ! lecture conditions limites
527            CALL interfoce_lim(itime, dtime, jour, &            CALL interfoce_lim(itime, dtime, jour, &
528                 &  klon, nisurf, knon, knindex, &                 klon, nisurf, knon, knindex, &
529                 &  debut, &                 debut, &
530                 &  tsurf_new, pctsrf_new)                 tsurf_new, pctsrf_new)
531    
532            !IM cf LF            !IM cf LF
533            DO ii = 1, knon            DO ii = 1, knon
534               tsurf_new(ii) = tsurf(ii)               tsurf_new(ii) = tsurf(ii)
535               !IMbad IF (pctsrf_new(ii,nisurf) < EPSFRA) then               !IMbad IF (pctsrf_new(ii, nisurf) < EPSFRA) then
536               IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then               IF (pctsrf_new(knindex(ii), nisurf) < EPSFRA) then
537                  snow(ii) = 0.0                  snow(ii) = 0.0
538                  !IM cf LF/JLD         tsurf(ii) = RTT - 1.8                  !IM cf LF/JLD tsurf(ii) = RTT - 1.8
539                  tsurf_new(ii) = RTT - 1.8                  tsurf_new(ii) = RTT - 1.8
540                  IF (soil_model) tsoil(ii,:) = RTT -1.8                  IF (soil_model) tsoil(ii, :) = RTT -1.8
541               endif               endif
542            enddo            enddo
543    
544            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
545    
546            IF (soil_model) THEN            IF (soil_model) THEN
547               !IM cf LF/JLD        CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)               !IM cf LF/JLD CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux)
548               CALL soil(dtime, nisurf, knon,snow, tsurf_new, tsoil,soilcap, soilflux)               CALL soil(dtime, nisurf, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
549               cal(1:knon) = RCPD / soilcap(1:knon)               cal(1:knon) = RCPD / soilcap(1:knon)
550               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)               radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
551               dif_grnd = 0.               dif_grnd = 0.
552            ELSE            ELSE
553               dif_grnd = 1.0 / tau_gl               dif_grnd = 1.0 / tau_gl
# Line 679  CONTAINS Line 560  CONTAINS
560         ENDIF         ENDIF
561    
562         CALL calcul_fluxs( klon, knon, nisurf, dtime, &         CALL calcul_fluxs( klon, knon, nisurf, dtime, &
563              &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &              tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
564              &   precip_rain, precip_snow, snow, qsurf,  &              precip_rain, precip_snow, snow, qsurf, &
565              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
566              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              petAcoef, peqAcoef, petBcoef, peqBcoef, &
567              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)              tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
568         !  
569         !IM: flux entre l'ocean et la glace de mer pour le "slab" ocean         !IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
570         DO i = 1, knon         DO i = 1, knon
571            flux_g(i) = 0.0            flux_g(i) = 0.0
572            !  
573            !IM: faire dependre le coefficient de conduction de la glace de mer            !IM: faire dependre le coefficient de conduction de la glace de mer
574            !    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.            ! de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
575            !    actuel correspond a 3m de glace de mer, cf. L.Li            ! actuel correspond a 3m de glace de mer, cf. L.Li
576            !  
577            !      IF(1.EQ.0) THEN            ! IF(1.EQ.0) THEN
578            !       IF(siceh(i).GT.0.) THEN            ! IF(siceh(i).GT.0.) THEN
579            !        new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)            ! new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)
580            !       ELSE            ! ELSE
581            !        new_dif_grnd(i) = 0.            ! new_dif_grnd(i) = 0.
582            !       ENDIF            ! ENDIF
583            !      ENDIF !(1.EQ.0) THEN            ! ENDIF !(1.EQ.0) THEN
584            !  
585            IF (cal(i).GT.1.0e-15) flux_g(i)=(tsurf_new(i)-t_grnd) &            IF (cal(i).GT.1.0e-15) flux_g(i)=(tsurf_new(i)-t_grnd) &
586                 &                          * dif_grnd(i) *RCPD/cal(i)                 * dif_grnd(i) *RCPD/cal(i)
587            !   &                          * new_dif_grnd(i) *RCPD/cal(i)            ! & * new_dif_grnd(i) *RCPD/cal(i)
588            tmp_flux_g(knindex(i))=flux_g(i)            tmp_flux_g(knindex(i))=flux_g(i)
589            tmp_radsol(knindex(i))=radsol(i)            tmp_radsol(knindex(i))=radsol(i)
590         ENDDO         ENDDO
591    
592         IF (ocean /= 'couple') THEN         CALL fonte_neige( klon, knon, nisurf, dtime, &
593            CALL fonte_neige( klon, knon, nisurf, dtime, &              tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
594                 &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &              precip_rain, precip_snow, snow, qsol, &
595                 &   precip_rain, precip_snow, snow, qsol,  &              radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
596                 &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              petAcoef, peqAcoef, petBcoef, peqBcoef, &
597                 &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
598                 &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &              fqcalving, ffonte, run_off_lic_0)
599                 &   fqcalving,ffonte, run_off_lic_0)  
600           ! calcul albedo
           !     calcul albedo  
   
           CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))    
           WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.  
           zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))  
           alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &  
                &                    0.6 * (1.0-zfra(1:knon))  
           !!      alb_new(1 : knon) = 0.6  
        ENDIF  
601    
602         fder_prev = fder             CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)
603           WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
604           zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon)/(snow(1:knon)+10.0)))
605           alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
606                0.6 * (1.0-zfra(1:knon))
607    
608           fder_prev = fder
609         fder = fder_prev + dflux_s + dflux_l         fder = fder_prev + dflux_s + dflux_l
610    
611         iloc = maxloc(fder(1:klon))         iloc = maxloc(fder(1:klon))
612         if (check.and.fder(iloc(1))> 0.) then         if (check.and.fder(iloc(1))> 0.) then
613            WRITE(*,*)'**** Debug fder ****'            WRITE(*, *)'**** Debug fder ****'
614            WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))            WRITE(*, *)'max fder(', iloc(1), ') = ', fder(iloc(1))
615            WRITE(*,*)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &            WRITE(*, *)'fder_prev, dflux_s, dflux_l', fder_prev(iloc(1)), &
616                 &                        dflux_s(iloc(1)), dflux_l(iloc(1))                 dflux_s(iloc(1)), dflux_l(iloc(1))
617         endif         endif
 !!$      where(fder.gt.0.)  
 !!$        fder = 0.  
 !!$      endwhere  
   
        !  
        ! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean  
        !  
        if (ocean == 'couple') then  
618    
           cumul =.true.  
619    
620            call interfoce_cpl(itime, dtime, cumul, &         ! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
                & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &  
                & ocean, npas, nexca, debut, lafin, &  
                & swdown, sollw, precip_rain, precip_snow, evap, tsurf, &  
                & fluxlat, fluxsens, fder, albedo, taux, tauy, &  
                & windsp, &  
                & zmasq, &  
                & tsurf_new, alb_new, &  
                & pctsrf_new)  
        endif  
621    
622         z0_new = 0.002         z0_new = 0.002
623         z0_new = SQRT(z0_new**2+rugoro**2)         z0_new = SQRT(z0_new**2+rugoro**2)
# Line 765  CONTAINS Line 625  CONTAINS
625    
626      else if (nisurf == is_lic) then      else if (nisurf == is_lic) then
627    
628         if (check) write(*,*)'glacier, nisurf = ',nisurf         if (check) write(*, *)'glacier, nisurf = ', nisurf
629    
630         if (.not. allocated(run_off_lic)) then         if (.not. allocated(run_off_lic)) then
631            allocate(run_off_lic(knon), stat = error)            allocate(run_off_lic(knon), stat = error)
632            if (error /= 0) then            if (error /= 0) then
633               abort_message='Pb allocation run_off_lic'               abort_message='Pb allocation run_off_lic'
634               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
635            endif            endif
636            run_off_lic = 0.            run_off_lic = 0.
637         endif         endif
638         !  
639         ! Surface "glacier continentaux" appel a l'interface avec le sol         ! Surface "glacier continentaux" appel a l'interface avec le sol
640         !  
641         IF (soil_model) THEN         IF (soil_model) THEN
642            CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil,soilcap, soilflux)            CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux)
643            cal(1:knon) = RCPD / soilcap(1:knon)            cal(1:knon) = RCPD / soilcap(1:knon)
644            radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)            radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
645         ELSE         ELSE
646            cal = RCPD * calice            cal = RCPD * calice
647            WHERE (snow > 0.0) cal = RCPD * calsno            WHERE (snow > 0.0) cal = RCPD * calsno
# Line 790  CONTAINS Line 650  CONTAINS
650         dif_grnd = 0.0         dif_grnd = 0.0
651    
652         call calcul_fluxs( klon, knon, nisurf, dtime, &         call calcul_fluxs( klon, knon, nisurf, dtime, &
653              &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &              tsurf, p1lay, cal, beta, tq_cdrag, ps, &
654              &   precip_rain, precip_snow, snow, qsurf,  &              precip_rain, precip_snow, snow, qsurf, &
655              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
656              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              petAcoef, peqAcoef, petBcoef, peqBcoef, &
657              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)              tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
658    
659         call fonte_neige( klon, knon, nisurf, dtime, &         call fonte_neige( klon, knon, nisurf, dtime, &
660              &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &              tsurf, p1lay, cal, beta, tq_cdrag, ps, &
661              &   precip_rain, precip_snow, snow, qsol,  &              precip_rain, precip_snow, snow, qsol, &
662              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
663              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              petAcoef, peqAcoef, petBcoef, peqBcoef, &
664              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &              tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
665              &   fqcalving,ffonte, run_off_lic_0)              fqcalving, ffonte, run_off_lic_0)
666    
667         ! passage du run-off des glaciers calcule dans fonte_neige au coupleur         ! passage du run-off des glaciers calcule dans fonte_neige au coupleur
668         bidule=0.         bidule=0.
669         bidule(1:knon)= run_off_lic(1:knon)             bidule(1:knon)= run_off_lic(1:knon)
670         call gath2cpl(bidule, tmp_rlic, klon, knon,iim,jjm,knindex)         call gath2cpl(bidule, tmp_rlic, klon, knon, iim, jjm, knindex)
671         !  
672         ! calcul albedo         ! calcul albedo
673         !  
674         CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))           CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)
675         WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.         WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
676         zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))         zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon)/(snow(1:knon)+10.0)))
677         alb_new(1 : knon)  = alb_neig(1 : knon)*zfra(1:knon) + &         alb_new(1 : knon) = alb_neig(1 : knon)*zfra(1:knon) + &
678              &                     0.6 * (1.0-zfra(1:knon))              0.6 * (1.0-zfra(1:knon))
679         !  
680         !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"         !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
681         !       alb_new(1 : knon)  = 0.6 !IM cf FH/GK         ! alb_new(1 : knon) = 0.6 !IM cf FH/GK
682         !       alb_new(1 : knon)  = 0.82         ! alb_new(1 : knon) = 0.82
683         !       alb_new(1 : knon)  = 0.77 !211003 Ksta0.77         ! alb_new(1 : knon) = 0.77 !211003 Ksta0.77
684         !       alb_new(1 : knon)  = 0.8 !KstaTER0.8 & LMD_ARMIP5         ! alb_new(1 : knon) = 0.8 !KstaTER0.8 & LMD_ARMIP5
685         !IM: KstaTER0.77 & LMD_ARMIP6             !IM: KstaTER0.77 & LMD_ARMIP6
686         alb_new(1 : knon)  = 0.77         alb_new(1 : knon) = 0.77
687    
688    
        !  
689         ! Rugosite         ! Rugosite
690         !  
691         z0_new = rugoro         z0_new = rugoro
692         !  
693         ! Remplissage des pourcentages de surface         ! Remplissage des pourcentages de surface
694         !  
695         pctsrf_new(:,nisurf) = pctsrf(:,nisurf)         pctsrf_new(:, nisurf) = pctsrf(:, nisurf)
696    
697         alblw(1:knon) = alb_new(1:knon)         alblw(1:knon) = alb_new(1:knon)
698      else      else
699         write(*,*)'Index surface = ',nisurf         write(*, *)'Index surface = ', nisurf
700         abort_message = 'Index surface non valable'         abort_message = 'Index surface non valable'
701         call abort_gcm(modname,abort_message,1)         call abort_gcm(modname, abort_message, 1)
702      endif      endif
703    
704    END SUBROUTINE interfsurf_hq    END SUBROUTINE interfsurf_hq
705    
706    !************************    !************************
707    
708    SUBROUTINE interfoce_cpl(itime, dtime, cumul, &    SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour,  &
709         & klon, iim, jjm, nisurf, pctsrf, knon, knindex, rlon, rlat, &         radsol, fluxo, fluxg, pctsrf,  &
710         & ocean, npas, nexca, debut, lafin, &         tslab, seaice, pctsrf_slab)
        & swdown, lwdown, precip_rain, precip_snow, evap, tsurf, &  
        & fluxlat, fluxsens, fder, albsol, taux, tauy, &  
        & windsp, &  
        & zmasq, &  
        & tsurf_new, alb_new, &  
        & pctsrf_new)  
   
     ! Cette routine sert d'interface entre le modele atmospherique et un  
     ! coupleur avec un modele d'ocean 'complet' derriere  
     !  
     ! Le modele de glace qu'il est prevu d'utiliser etant couple directement a  
     ! l'ocean presentement, on va passer deux fois dans cette routine par pas de  
     ! temps physique, une fois avec les points oceans et l'autre avec les points  
     ! glace. A chaque pas de temps de couplage, la lecture des champs provenant  
     ! du coupleur se fera "dans" l'ocean et l'ecriture des champs a envoyer  
     ! au coupleur "dans" la glace. Il faut donc des tableaux de travail "tampons"  
     ! dimensionnes sur toute la grille qui remplissent les champs sur les  
     ! domaines ocean/glace quand il le faut. Il est aussi necessaire que l'index  
     ! ocean soit traiter avant l'index glace (sinon tout intervertir)  
     !  
     !  
     ! L. Fairhead 02/2000  
     !  
     ! input:  
     !   itime        numero du pas de temps  
     !   iim, jjm     nbres de pts de grille  
     !   dtime        pas de temps de la physique  
     !   klon         nombre total de points de grille  
     !   nisurf       index de la surface a traiter (1 = sol continental)  
     !   pctsrf       tableau des fractions de surface de chaque maille  
     !   knon         nombre de points de la surface a traiter  
     !   knindex      index des points de la surface a traiter  
     !   rlon         longitudes  
     !   rlat         latitudes  
     !   debut        logical: 1er appel a la physique  
     !   lafin        logical: dernier appel a la physique  
     !   ocean        type d'ocean  
     !   nexca        frequence de couplage  
     !   swdown       flux solaire entrant a la surface  
     !   lwdown       flux IR net a la surface  
     !   precip_rain  precipitation liquide  
     !   precip_snow  precipitation solide  
     !   evap         evaporation  
     !   tsurf        temperature de surface  
     !   fder         derivee dF/dT  
     !   albsol       albedo du sol (coherent avec swdown)  
     !   taux         tension de vent en x  
     !   tauy         tension de vent en y  
     !    windsp       module du vent a 10m  
     !   nexca        frequence de couplage  
     !   zmasq        masque terre/ocean  
     !  
     !  
     ! output:  
     !   tsurf_new    temperature au sol  
     !   alb_new      albedo  
     !   pctsrf_new   nouvelle repartition des surfaces  
     !   alb_ice      albedo de la glace  
     !  
     use temps  
     use iniprint  
     use abort_gcm_m, only: abort_gcm  
     use gath_cpl, only: gath2cpl, cpl2gath  
     use ioipsl  
     use indicesol  
     use YOMCST  
   
     ! Parametres d'entree  
     integer, intent(IN) :: itime  
     integer, intent(IN) :: iim, jjm  
     real, intent(IN) :: dtime  
     integer, intent(IN) :: klon  
     integer, intent(IN) :: nisurf  
     integer, intent(IN) :: knon  
     real, dimension(klon,nbsrf), intent(IN) :: pctsrf  
     integer, dimension(klon), intent(in) :: knindex  
     logical, intent(IN) :: debut, lafin  
     real, dimension(klon), intent(IN) :: rlon, rlat  
     character (len = 6)  :: ocean  
     real, dimension(klon), intent(IN) :: lwdown, swdown  
     real, dimension(klon), intent(IN) :: precip_rain, precip_snow  
     real, dimension(klon), intent(IN) :: tsurf, fder, albsol, taux, tauy  
     real, dimension(klon), intent(IN) :: windsp  
     INTEGER              :: nexca, npas  
     real, dimension(klon), intent(IN) :: zmasq  
     real, dimension(klon), intent(IN) :: fluxlat, fluxsens  
     logical, intent(IN)               :: cumul  
     real, dimension(klon), intent(INOUT) :: evap  
   
     ! Parametres de sortie  
     real, dimension(klon), intent(OUT):: tsurf_new, alb_new  
     real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new  
   
     ! Variables locales  
     integer                    :: j, error, sum_error, ig, cpl_index,i  
     INTEGER :: nsrf  
     character (len = 20) :: modname = 'interfoce_cpl'  
     character (len = 80) :: abort_message  
     logical,save              :: check = .FALSE.  
     ! variables pour moyenner les variables de couplage  
     real, allocatable, dimension(:,:),save :: cpl_sols, cpl_nsol, cpl_rain  
     real, allocatable, dimension(:,:),save :: cpl_snow, cpl_evap, cpl_tsol  
     real, allocatable, dimension(:,:),save :: cpl_fder, cpl_albe, cpl_taux  
     real, allocatable, dimension(:,:),save :: cpl_windsp  
     real, allocatable, dimension(:,:),save :: cpl_tauy  
     REAL, ALLOCATABLE, DIMENSION(:,:),SAVE :: cpl_rriv, cpl_rcoa, cpl_rlic  
 !!$  
     ! variables tampons avant le passage au coupleur  
     real, allocatable, dimension(:,:,:),save :: tmp_sols, tmp_nsol, tmp_rain  
     real, allocatable, dimension(:,:,:),save :: tmp_snow, tmp_evap, tmp_tsol  
     real, allocatable, dimension(:,:,:),save :: tmp_fder, tmp_albe, tmp_taux  
     real, allocatable, dimension(:,:,:),save :: tmp_windsp  
     REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: tmp_tauy  
     ! variables a passer au coupleur  
     real, dimension(iim, jjm+1) :: wri_sol_ice, wri_sol_sea, wri_nsol_ice  
     real, dimension(iim, jjm+1) :: wri_nsol_sea, wri_fder_ice, wri_evap_ice  
     REAL, DIMENSION(iim, jjm+1) :: wri_evap_sea, wri_rcoa, wri_rriv  
     REAL, DIMENSION(iim, jjm+1) :: wri_rain, wri_snow, wri_taux, wri_tauy  
     REAL, DIMENSION(iim, jjm+1) :: wri_windsp  
     REAL, DIMENSION(iim, jjm+1) :: wri_calv  
     REAL, DIMENSION(iim, jjm+1) :: wri_tauxx, wri_tauyy, wri_tauzz  
     REAL, DIMENSION(iim, jjm+1) :: tmp_lon, tmp_lat  
     ! variables relues par le coupleur  
     ! read_sic = fraction de glace  
     ! read_sit = temperature de glace  
     real, allocatable, dimension(:,:),save :: read_sst, read_sic, read_sit  
     real, allocatable, dimension(:,:),save :: read_alb_sic  
     ! variable tampon  
     real, dimension(klon)       :: tamp_sic  
     ! sauvegarde des fractions de surface d'un pas de temps a l'autre apres  
     ! l'avoir lu  
     real, allocatable,dimension(:,:),save :: pctsrf_sav  
     real, dimension(iim, jjm+1, 2) :: tamp_srf  
     integer, allocatable, dimension(:), save :: tamp_ind  
     real, allocatable, dimension(:,:),save :: tamp_zmasq  
     real, dimension(iim, jjm+1) :: deno  
     integer                     :: idtime  
     integer, allocatable,dimension(:),save :: unity  
     !  
     logical, save    :: first_appel = .true.  
     logical,save          :: print  
     !maf  
     ! variables pour avoir une sortie IOIPSL des champs echanges  
     CHARACTER(len=80),SAVE :: clintocplnam, clfromcplnam  
     INTEGER, SAVE :: jf,nhoridct,nidct  
     INTEGER, SAVE :: nhoridcs,nidcs  
     INTEGER :: ndexct(iim*(jjm+1)),ndexcs(iim*(jjm+1))  
     REAL :: zx_lon(iim,jjm+1), zx_lat(iim,jjm+1), zjulian  
     INTEGER,save :: idayref  
     !med  integer :: itau_w  
     integer,save :: itau_w  
     integer :: nb_interf_cpl  
     include "param_cou.h"  
     include "inc_cpl.h"  
     !  
     ! Initialisation  
     !  
     if (check) write(*,*)'Entree ',modname,'nisurf = ',nisurf  
   
     if (first_appel) then  
        error = 0  
        allocate(unity(klon), stat = error)  
        if ( error  /=0) then  
           abort_message='Pb allocation variable unity'  
           call abort_gcm(modname,abort_message,1)  
        endif  
        allocate(pctsrf_sav(klon,nbsrf), stat = error)  
        if ( error  /=0) then  
           abort_message='Pb allocation variable pctsrf_sav'  
           call abort_gcm(modname,abort_message,1)  
        endif  
        pctsrf_sav = 0.  
   
        do ig = 1, klon  
           unity(ig) = ig  
        enddo  
        sum_error = 0  
        allocate(cpl_sols(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_nsol(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_rain(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_snow(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_evap(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_tsol(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_fder(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_albe(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_taux(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_windsp(klon,2), stat = error); sum_error = sum_error + error  
        allocate(cpl_tauy(klon,2), stat = error); sum_error = sum_error + error  
        ALLOCATE(cpl_rriv(iim,jjm+1), stat=error); sum_error = sum_error + error  
        ALLOCATE(cpl_rcoa(iim,jjm+1), stat=error); sum_error = sum_error + error  
        ALLOCATE(cpl_rlic(iim,jjm+1), stat=error); sum_error = sum_error + error  
        !!  
        allocate(read_sst(iim, jjm+1), stat = error); sum_error = sum_error + error  
        allocate(read_sic(iim, jjm+1), stat = error); sum_error = sum_error + error  
        allocate(read_sit(iim, jjm+1), stat = error); sum_error = sum_error + error  
        allocate(read_alb_sic(iim, jjm+1), stat = error); sum_error = sum_error + error  
   
        if (sum_error /= 0) then  
           abort_message='Pb allocation variables couplees'  
           call abort_gcm(modname,abort_message,1)  
        endif  
        cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.  
        cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.  
        cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.; cpl_rlic = 0.  
        cpl_windsp = 0.  
   
        sum_error = 0  
        allocate(tamp_ind(klon), stat = error); sum_error = sum_error + error  
        allocate(tamp_zmasq(iim, jjm+1), stat = error); sum_error = sum_error + error      
        do ig = 1, klon  
           tamp_ind(ig) = ig  
        enddo  
        call gath2cpl(zmasq, tamp_zmasq, klon, klon, iim, jjm, tamp_ind)  
        !  
        ! initialisation couplage  
        !  
        idtime = int(dtime)  
        !  
        ! initialisation sorties netcdf  
        !  
        idayref = day_ini  
        CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)  
        CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,zx_lon)  
        DO i = 1, iim  
           zx_lon(i,1) = rlon(i+1)  
           zx_lon(i,jjm+1) = rlon(i+1)  
        ENDDO  
        CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,zx_lat)  
        clintocplnam="cpl_atm_tauflx"  
        CALL histbeg_totreg(clintocplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &  
             & itau_phy,zjulian,dtime,nhoridct,nidct)  
        ! no vertical axis  
        CALL histdef(nidct, 'tauxe','tauxe', &  
             & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)  
        CALL histdef(nidct, 'tauyn','tauyn', &  
             & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)  
        CALL histdef(nidct, 'tmp_lon','tmp_lon', &  
             & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)  
        CALL histdef(nidct, 'tmp_lat','tmp_lat', &  
             & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)  
        DO jf=1,jpflda2o1 + jpflda2o2  
           CALL histdef(nidct, cl_writ(jf),cl_writ(jf), &  
                & "-",iim, jjm+1, nhoridct, 1, 1, 1, -99, 32, "inst", dtime,dtime)  
        END DO  
        CALL histend(nidct)  
        CALL histsync(nidct)  
   
        clfromcplnam="cpl_atm_sst"  
        CALL histbeg_totreg(clfromcplnam, iim,zx_lon(:,1),jjm+1,zx_lat(1,:),1,iim,1,jjm+1, &  
             & 0,zjulian,dtime,nhoridcs,nidcs)  
        ! no vertical axis  
        DO jf=1,jpfldo2a  
           CALL histdef(nidcs, cl_read(jf),cl_read(jf), &  
                & "-",iim, jjm+1, nhoridcs, 1, 1, 1, -99, 32, "inst", dtime,dtime)  
        END DO  
        CALL histend(nidcs)  
        CALL histsync(nidcs)  
   
        ! pour simuler la fonte des glaciers antarctiques  
        !  
        surf_maille = (4. * rpi * ra**2) / (iim * (jjm +1))  
        ALLOCATE(coeff_iceberg(iim,jjm+1), stat=error)  
        if (error /= 0) then  
           abort_message='Pb allocation variable coeff_iceberg'  
           call abort_gcm(modname,abort_message,1)  
        endif  
        open (12,file='flux_iceberg',form='formatted',status='old')  
        read (12,*) coeff_iceberg  
        close (12)  
        num_antarctic = max(1, count(coeff_iceberg > 0))  
   
        first_appel = .false.  
     endif ! fin if (first_appel)  
   
     ! Initialisations  
   
     ! calcul des fluxs a passer  
     nb_interf_cpl = nb_interf_cpl + 1    
     if (check) write(lunout,*)'passage dans interface_surf.F90 :  ',nb_interf_cpl  
     cpl_index = 1  
     if (nisurf == is_sic) cpl_index = 2  
     if (cumul) then  
        if (check) write(lunout,*)'passage dans cumul '  
        if (check) write(lunout,*)'valeur de cpl_index ', cpl_index  
        ! -- LOOP    
        if (check) write(*,*) modname, 'cumul des champs'  
        do ig = 1, knon  
           cpl_sols(ig,cpl_index) = cpl_sols(ig,cpl_index) &  
                &                          + swdown(ig)      / FLOAT(nexca)  
           cpl_nsol(ig,cpl_index) = cpl_nsol(ig,cpl_index) &  
                &                          + (lwdown(ig) + fluxlat(ig) +fluxsens(ig))&  
                &                                / FLOAT(nexca)  
           cpl_rain(ig,cpl_index) = cpl_rain(ig,cpl_index) &  
                &                          + precip_rain(ig) / FLOAT(nexca)  
           cpl_snow(ig,cpl_index) = cpl_snow(ig,cpl_index) &  
                &                          + precip_snow(ig) / FLOAT(nexca)  
           cpl_evap(ig,cpl_index) = cpl_evap(ig,cpl_index) &  
                &                          + evap(ig)        / FLOAT(nexca)  
           cpl_tsol(ig,cpl_index) = cpl_tsol(ig,cpl_index) &  
                &                          + tsurf(ig)       / FLOAT(nexca)  
           cpl_fder(ig,cpl_index) = cpl_fder(ig,cpl_index) &  
                &                          + fder(ig)        / FLOAT(nexca)  
           cpl_albe(ig,cpl_index) = cpl_albe(ig,cpl_index) &  
                &                          + albsol(ig)      / FLOAT(nexca)  
           cpl_taux(ig,cpl_index) = cpl_taux(ig,cpl_index) &  
                &                          + taux(ig)        / FLOAT(nexca)  
           cpl_tauy(ig,cpl_index) = cpl_tauy(ig,cpl_index) &  
                &                          + tauy(ig)        / FLOAT(nexca)  
           IF (cpl_index .EQ. 1) THEN  
              cpl_windsp(ig,cpl_index) = cpl_windsp(ig,cpl_index) &  
                   &                          + windsp(ig)      / FLOAT(nexca)  
           ENDIF  
        enddo  
        IF (cpl_index .EQ. 1) THEN  
           cpl_rriv(:,:) = cpl_rriv(:,:) + tmp_rriv(:,:) / FLOAT(nexca)  
           cpl_rcoa(:,:) = cpl_rcoa(:,:) + tmp_rcoa(:,:) / FLOAT(nexca)  
           cpl_rlic(:,:) = cpl_rlic(:,:) + tmp_rlic(:,:) / FLOAT(nexca)  
        ENDIF  
     endif  
   
     if (mod(itime, nexca) == 1) then  
        !  
        ! Demande des champs au coupleur  
        !  
        ! Si le domaine considere est l'ocean, on lit les champs venant du coupleur  
        !  
        if (nisurf == is_oce .and. .not. cumul) then  
           if (check) write(*,*)'rentree fromcpl, itime-1 = ',itime-1  
           !  
           ! sorties NETCDF des champs recus  
           !  
           ndexcs(:)=0  
           itau_w = itau_phy + itime  
           CALL histwrite(nidcs,cl_read(1),itau_w,read_sst,iim*(jjm+1),ndexcs)  
           CALL histwrite(nidcs,cl_read(2),itau_w,read_sic,iim*(jjm+1),ndexcs)  
           CALL histwrite(nidcs,cl_read(3),itau_w,read_alb_sic,iim*(jjm+1),ndexcs)  
           CALL histwrite(nidcs,cl_read(4),itau_w,read_sit,iim*(jjm+1),ndexcs)  
           CALL histsync(nidcs)  
           ! pas utile      IF (npas-itime.LT.nexca )CALL histclo(nidcs)  
   
           do j = 1, jjm + 1  
              do ig = 1, iim  
                 if (abs(1. - read_sic(ig,j)) < 0.00001) then  
                    read_sst(ig,j) = RTT - 1.8  
                    read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)  
                    read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)  
                 else if (abs(read_sic(ig,j)) < 0.00001) then  
                    read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))  
                    read_sit(ig,j) = read_sst(ig,j)  
                    read_alb_sic(ig,j) =  0.6  
                 else  
                    read_sst(ig,j) = read_sst(ig,j) / (1. - read_sic(ig,j))  
                    read_sit(ig,j) = read_sit(ig,j) / read_sic(ig,j)  
                    read_alb_sic(ig,j) = read_alb_sic(ig,j) / read_sic(ig,j)  
                 endif  
              enddo  
           enddo  
           !  
           ! transformer read_sic en pctsrf_sav  
           !  
           call cpl2gath(read_sic, tamp_sic , klon, klon,iim,jjm, unity)  
           do ig = 1, klon  
              IF (pctsrf(ig,is_oce) > epsfra .OR.            &  
                   &             pctsrf(ig,is_sic) > epsfra) THEN  
                 pctsrf_sav(ig,is_sic) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &  
                      &                               * tamp_sic(ig)  
                 pctsrf_sav(ig,is_oce) = (pctsrf(ig,is_oce) + pctsrf(ig,is_sic)) &  
                      &                        - pctsrf_sav(ig,is_sic)  
              endif  
           enddo  
           !  
           ! Pour rattraper des erreurs d'arrondis  
           !  
           where (abs(pctsrf_sav(:,is_sic)) .le. 2.*epsilon(pctsrf_sav(1,is_sic)))  
              pctsrf_sav(:,is_sic) = 0.  
              pctsrf_sav(:,is_oce) = pctsrf(:,is_oce) + pctsrf(:,is_sic)  
           endwhere  
           where (abs(pctsrf_sav(:,is_oce)) .le. 2.*epsilon(pctsrf_sav(1,is_oce)))  
              pctsrf_sav(:,is_sic) = pctsrf(:,is_oce) + pctsrf(:,is_sic)  
              pctsrf_sav(:,is_oce) = 0.  
           endwhere  
           if (minval(pctsrf_sav(:,is_oce)) < 0.) then  
              write(*,*)'Pb fraction ocean inferieure a 0'  
              write(*,*)'au point ',minloc(pctsrf_sav(:,is_oce))  
              write(*,*)'valeur = ',minval(pctsrf_sav(:,is_oce))  
              abort_message = 'voir ci-dessus'  
              call abort_gcm(modname,abort_message,1)  
           endif  
           if (minval(pctsrf_sav(:,is_sic)) < 0.) then  
              write(*,*)'Pb fraction glace inferieure a 0'  
              write(*,*)'au point ',minloc(pctsrf_sav(:,is_sic))  
              write(*,*)'valeur = ',minval(pctsrf_sav(:,is_sic))  
              abort_message = 'voir ci-dessus'  
              call abort_gcm(modname,abort_message,1)  
           endif  
        endif  
     endif                         ! fin mod(itime, nexca) == 1  
   
     if (mod(itime, nexca) == 0) then  
        !  
        ! allocation memoire  
        if (nisurf == is_oce .and. (.not. cumul) ) then  
           sum_error = 0  
           allocate(tmp_sols(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_nsol(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_rain(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_snow(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_evap(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_tsol(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_fder(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_albe(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_taux(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_tauy(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           allocate(tmp_windsp(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
 !!$      allocate(tmp_rriv(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
 !!$      allocate(tmp_rcoa(iim,jjm+1,2), stat=error); sum_error = sum_error + error  
           if (sum_error /= 0) then  
              abort_message='Pb allocation variables couplees pour l''ecriture'  
              call abort_gcm(modname,abort_message,1)  
           endif  
        endif  
   
        !  
        ! Mise sur la bonne grille des champs a passer au coupleur  
        !  
        cpl_index = 1  
        if (nisurf == is_sic) cpl_index = 2  
        call gath2cpl(cpl_sols(1,cpl_index), tmp_sols(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_nsol(1,cpl_index), tmp_nsol(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_rain(1,cpl_index), tmp_rain(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_snow(1,cpl_index), tmp_snow(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_evap(1,cpl_index), tmp_evap(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_tsol(1,cpl_index), tmp_tsol(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_fder(1,cpl_index), tmp_fder(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_albe(1,cpl_index), tmp_albe(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_taux(1,cpl_index), tmp_taux(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
        call gath2cpl(cpl_windsp(1,cpl_index), tmp_windsp(1,1,cpl_index), klon, knon,iim,jjm,             knindex)  
        call gath2cpl(cpl_tauy(1,cpl_index), tmp_tauy(1,1,cpl_index), klon, knon,iim,jjm,                  knindex)  
   
        !  
        ! Si le domaine considere est la banquise, on envoie les champs au coupleur  
        !  
        if (nisurf == is_sic .and. cumul) then  
           wri_rain = 0.; wri_snow = 0.; wri_rcoa = 0.; wri_rriv = 0.  
           wri_taux = 0.; wri_tauy = 0.  
           wri_windsp = 0.  
           ! -- LOOP        
           call gath2cpl(pctsrf(1,is_oce), tamp_srf(1,1,1), klon, klon, iim, jjm, tamp_ind)  
           call gath2cpl(pctsrf(1,is_sic), tamp_srf(1,1,2), klon, klon, iim, jjm, tamp_ind)  
   
           wri_sol_ice = tmp_sols(:,:,2)  
           wri_sol_sea = tmp_sols(:,:,1)  
           wri_nsol_ice = tmp_nsol(:,:,2)  
           wri_nsol_sea = tmp_nsol(:,:,1)  
           wri_fder_ice = tmp_fder(:,:,2)  
           wri_evap_ice = tmp_evap(:,:,2)  
           wri_evap_sea = tmp_evap(:,:,1)  
           wri_windsp = tmp_windsp(:,:,1)  
711    
 !!$PB  
           wri_rriv = cpl_rriv(:,:)  
           wri_rcoa = cpl_rcoa(:,:)  
           DO j = 1, jjm + 1  
              wri_calv(:,j) = sum(cpl_rlic(:,j)) / iim  
           enddo  
   
           where (tamp_zmasq /= 1.)  
              deno =  tamp_srf(:,:,1) + tamp_srf(:,:,2)  
              wri_rain = tmp_rain(:,:,1) * tamp_srf(:,:,1) / deno +    &  
                   &            tmp_rain(:,:,2) * tamp_srf(:,:,2) / deno  
              wri_snow = tmp_snow(:,:,1) * tamp_srf(:,:,1) / deno +    &  
                   &            tmp_snow(:,:,2) * tamp_srf(:,:,2) / deno  
              wri_taux = tmp_taux(:,:,1) * tamp_srf(:,:,1) / deno +    &  
                   &            tmp_taux(:,:,2) * tamp_srf(:,:,2) / deno  
              wri_tauy = tmp_tauy(:,:,1) * tamp_srf(:,:,1) / deno +    &  
                   &            tmp_tauy(:,:,2) * tamp_srf(:,:,2) / deno  
           endwhere  
           !  
           ! pour simuler la fonte des glaciers antarctiques  
           !  
           !$$$        wri_rain = wri_rain      &  
           !$$$      &     + coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)  
           !      wri_calv = coeff_iceberg * cte_flux_iceberg / (num_antarctic * surf_maille)  
           !  
           ! on passe les coordonnées de la grille  
           !  
   
           CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlon,tmp_lon)  
           CALL gr_fi_ecrit(1,klon,iim,jjm+1,rlat,tmp_lat)  
   
           DO i = 1, iim  
              tmp_lon(i,1) = rlon(i+1)  
              tmp_lon(i,jjm + 1) = rlon(i+1)  
           ENDDO  
           !  
           ! sortie netcdf des champs pour le changement de repere  
           !  
           ndexct(:)=0  
           CALL histwrite(nidct,'tauxe',itau_w,wri_taux,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,'tauyn',itau_w,wri_tauy,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,'tmp_lon',itau_w,tmp_lon,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,'tmp_lat',itau_w,tmp_lat,iim*(jjm+1),ndexct)  
   
           !  
           ! calcul 3 coordonnées du vent  
           !  
           CALL atm2geo (iim , jjm + 1, wri_taux, wri_tauy, tmp_lon, tmp_lat, &  
                & wri_tauxx, wri_tauyy, wri_tauzz )  
           !  
           ! sortie netcdf des champs apres changement de repere et juste avant  
           ! envoi au coupleur  
           !  
           CALL histwrite(nidct,cl_writ(8),itau_w,wri_sol_ice,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(9),itau_w,wri_sol_sea,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(10),itau_w,wri_nsol_ice,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(11),itau_w,wri_nsol_sea,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(12),itau_w,wri_fder_ice,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(13),itau_w,wri_evap_ice,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(14),itau_w,wri_evap_sea,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(15),itau_w,wri_rain,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(16),itau_w,wri_snow,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(17),itau_w,wri_rcoa,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(18),itau_w,wri_rriv,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(19),itau_w,wri_calv,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(1),itau_w,wri_tauxx,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(2),itau_w,wri_tauyy,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(3),itau_w,wri_tauzz,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(4),itau_w,wri_tauxx,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(5),itau_w,wri_tauyy,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(6),itau_w,wri_tauzz,iim*(jjm+1),ndexct)  
           CALL histwrite(nidct,cl_writ(7),itau_w,wri_windsp,iim*(jjm+1),ndexct)  
           CALL histsync(nidct)  
           ! pas utile      IF (lafin) CALL histclo(nidct)  
           !  
           cpl_sols = 0.; cpl_nsol = 0.; cpl_rain = 0.; cpl_snow = 0.  
           cpl_evap = 0.; cpl_tsol = 0.; cpl_fder = 0.; cpl_albe = 0.  
           cpl_taux = 0.; cpl_tauy = 0.; cpl_rriv = 0.; cpl_rcoa = 0.; cpl_rlic = 0.  
           cpl_windsp = 0.  
           !  
           ! deallocation memoire variables temporaires  
           !  
           sum_error = 0  
           deallocate(tmp_sols, stat=error); sum_error = sum_error + error  
           deallocate(tmp_nsol, stat=error); sum_error = sum_error + error  
           deallocate(tmp_rain, stat=error); sum_error = sum_error + error  
           deallocate(tmp_snow, stat=error); sum_error = sum_error + error  
           deallocate(tmp_evap, stat=error); sum_error = sum_error + error  
           deallocate(tmp_fder, stat=error); sum_error = sum_error + error  
           deallocate(tmp_tsol, stat=error); sum_error = sum_error + error  
           deallocate(tmp_albe, stat=error); sum_error = sum_error + error  
           deallocate(tmp_taux, stat=error); sum_error = sum_error + error  
           deallocate(tmp_tauy, stat=error); sum_error = sum_error + error  
           deallocate(tmp_windsp, stat=error); sum_error = sum_error + error  
           if (sum_error /= 0) then  
              abort_message='Pb deallocation variables couplees'  
              call abort_gcm(modname,abort_message,1)  
           endif  
   
        endif  
   
     endif            ! fin (mod(itime, nexca) == 0)  
     !  
     ! on range les variables lues/sauvegardees dans les bonnes variables de sortie  
     !  
     if (nisurf == is_oce) then  
        call cpl2gath(read_sst, tsurf_new, klon, knon,iim,jjm, knindex)  
     else if (nisurf == is_sic) then  
        call cpl2gath(read_sit, tsurf_new, klon, knon,iim,jjm, knindex)  
        call cpl2gath(read_alb_sic, alb_new, klon, knon,iim,jjm, knindex)  
     endif  
     pctsrf_new(:,nisurf) = pctsrf_sav(:,nisurf)  
   
     !  if (lafin) call quitcpl  
   
   END SUBROUTINE interfoce_cpl  
   
   !************************  
   
   SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &  
        & radsol, fluxo, fluxg, pctsrf, &  
        & tslab, seaice, pctsrf_slab)  
     !  
712      ! Cette routine calcule la temperature d'un slab ocean, la glace de mer      ! Cette routine calcule la temperature d'un slab ocean, la glace de mer
713      ! et les pourcentages de la maille couverte par l'ocean libre et/ou      ! et les pourcentages de la maille couverte par l'ocean libre et/ou
714      ! la glace de mer pour un "slab" ocean de 50m      ! la glace de mer pour un "slab" ocean de 50m
715      !  
716      ! I. Musat 04.02.2005      ! I. Musat 04.02.2005
717      !  
718      ! input:      ! input:
719      !   klon         nombre total de points de grille      ! klon nombre total de points de grille
720      !   debut        logical: 1er appel a la physique      ! debut logical: 1er appel a la physique
721      !   itap         numero du pas de temps      ! itap numero du pas de temps
722      !   dtime        pas de temps de la physique (en s)      ! dtime pas de temps de la physique (en s)
723      !   ijour        jour dans l'annee en cours      ! ijour jour dans l'annee en cours
724      !   radsol       rayonnement net au sol (LW + SW)      ! radsol rayonnement net au sol (LW + SW)
725      !   fluxo        flux turbulent (sensible + latent) sur les mailles oceaniques      ! fluxo flux turbulent (sensible + latent) sur les mailles oceaniques
726      !   fluxg        flux de conduction entre la surface de la glace de mer et l'ocean      ! fluxg flux de conduction entre la surface de la glace de mer et l'ocean
727      !   pctsrf       tableau des pourcentages de surface de chaque maille      ! pctsrf tableau des pourcentages de surface de chaque maille
728      ! output:      ! output:
729      !   tslab        temperature de l'ocean libre      ! tslab temperature de l'ocean libre
730      !   seaice       glace de mer (kg/m2)      ! seaice glace de mer (kg/m2)
731      !   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab      ! pctsrf_slab "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
732      !  
733      use indicesol      use indicesol
734      use clesphys      use clesphys
735      use abort_gcm_m, only: abort_gcm      use abort_gcm_m, only: abort_gcm
736      use YOMCST      use SUPHEC_M
737    
738      ! Parametres d'entree      ! Parametres d'entree
739      integer, intent(IN) :: klon      integer, intent(IN) :: klon
# Line 1468  CONTAINS Line 747  CONTAINS
747      real, dimension(klon, nbsrf), intent(IN) :: pctsrf      real, dimension(klon, nbsrf), intent(IN) :: pctsrf
748      ! Parametres de sortie      ! Parametres de sortie
749      real, dimension(klon), intent(INOUT) :: tslab      real, dimension(klon), intent(INOUT) :: tslab
750      real, dimension(klon), intent(INOUT)        :: seaice ! glace de mer (kg/m2)      real, dimension(klon), intent(INOUT) :: seaice ! glace de mer (kg/m2)
751      real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_slab      real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_slab
752      !  
753      ! Variables locales :      ! Variables locales :
754      INTEGER, save :: lmt_pas, julien, idayvrai      INTEGER, save :: lmt_pas, julien, idayvrai
755      REAL, parameter :: unjour=86400.      REAL, parameter :: unjour=86400.
756      real, allocatable, dimension(:), save :: tmp_tslab, tmp_seaice      real, allocatable, dimension(:), save :: tmp_tslab, tmp_seaice
757      REAL, allocatable, dimension(:), save :: slab_bils      REAL, allocatable, dimension(:), save :: slab_bils
758      REAL, allocatable, dimension(:), save :: lmt_bils      REAL, allocatable, dimension(:), save :: lmt_bils
759      logical,save              :: check = .false.      logical, save :: check = .false.
760      !  
761      REAL, parameter :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)      REAL, parameter :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
762      REAL, parameter :: cbing=0.334e+05        ! J/kg      REAL, parameter :: cbing=0.334e+05 ! J/kg
763      real, dimension(klon)                 :: siceh !hauteur de la glace de mer (m)      real, dimension(klon) :: siceh !hauteur de la glace de mer (m)
764      INTEGER :: i      INTEGER :: i
765      integer :: sum_error, error      integer :: sum_error, error
766      REAL :: zz, za, zb      REAL :: zz, za, zb
767      !  
768      character (len = 80) :: abort_message      character (len = 80) :: abort_message
769      character (len = 20) :: modname = 'interfoce_slab'      character (len = 20) :: modname = 'interfoce_slab'
770      !  
771      julien = MOD(ijour,360)      julien = MOD(ijour, 360)
772      sum_error = 0      sum_error = 0
773      IF (debut) THEN      IF (debut) THEN
774         allocate(slab_bils(klon), stat = error); sum_error = sum_error + error         allocate(slab_bils(klon), stat = error); sum_error = sum_error + error
# Line 1497  CONTAINS Line 776  CONTAINS
776         allocate(tmp_tslab(klon), stat = error); sum_error = sum_error + error         allocate(tmp_tslab(klon), stat = error); sum_error = sum_error + error
777         allocate(tmp_seaice(klon), stat = error); sum_error = sum_error + error         allocate(tmp_seaice(klon), stat = error); sum_error = sum_error + error
778         if (sum_error /= 0) then         if (sum_error /= 0) then
779            abort_message='Pb allocation var. slab_bils,lmt_bils,tmp_tslab,tmp_seaice'            abort_message='Pb allocation var. slab_bils, lmt_bils, tmp_tslab, tmp_seaice'
780            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
781         endif         endif
782         tmp_tslab=tslab         tmp_tslab=tslab
783         tmp_seaice=seaice         tmp_seaice=seaice
784         lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour         lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
785         !  
786         IF (check) THEN         IF (check) THEN
787            PRINT*,'interfoce_slab klon, debut, itap, dtime, ijour, &            PRINT*, 'interfoce_slab klon, debut, itap, dtime, ijour, &
788                 &          lmt_pas ', klon, debut, itap, dtime, ijour, &                 & lmt_pas ', klon, debut, itap, dtime, ijour, &
789                 &          lmt_pas                 lmt_pas
790         ENDIF !check         ENDIF !check
791         !  
792         PRINT*, '************************'         PRINT*, '************************'
793         PRINT*, 'SLAB OCEAN est actif, prenez precautions !'         PRINT*, 'SLAB OCEAN est actif, prenez precautions !'
794         PRINT*, '************************'         PRINT*, '************************'
795         !  
796         ! a mettre un slab_bils aussi en force !!!         ! a mettre un slab_bils aussi en force !!!
797         !  
798         DO i = 1, klon         DO i = 1, klon
799            slab_bils(i) = 0.0            slab_bils(i) = 0.0
800         ENDDO         ENDDO
801         !  
802      ENDIF !debut      ENDIF !debut
803      pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)      pctsrf_slab(1:klon, 1:nbsrf) = pctsrf(1:klon, 1:nbsrf)
804      !  
805      ! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee      ! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
806      !  
807      IF (MOD(itap,lmt_pas) .EQ. 1) THEN !1er pas de temps de la journee      IF (MOD(itap, lmt_pas) .EQ. 1) THEN !1er pas de temps de la journee
808         idayvrai = ijour         idayvrai = ijour
809         CALL condsurf(julien,idayvrai, lmt_bils)         CALL condsurf(julien, idayvrai, lmt_bils)
810      ENDIF !(MOD(itap-1,lmt_pas) .EQ. 0) THEN      ENDIF !(MOD(itap-1, lmt_pas) .EQ. 0) THEN
811    
812      DO i = 1, klon      DO i = 1, klon
813         IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &         IF((pctsrf_slab(i, is_oce).GT.epsfra).OR. &
814              &  (pctsrf_slab(i,is_sic).GT.epsfra)) THEN              (pctsrf_slab(i, is_sic).GT.epsfra)) THEN
815            !  
816            ! fabriquer de la glace si congelation atteinte:            ! fabriquer de la glace si congelation atteinte:
817            !  
818            IF (tmp_tslab(i).LT.(RTT-1.8)) THEN            IF (tmp_tslab(i).LT.(RTT-1.8)) THEN
819               zz =  (RTT-1.8)-tmp_tslab(i)               zz = (RTT-1.8)-tmp_tslab(i)
820               tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz               tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz
821               seaice(i) = tmp_seaice(i)               seaice(i) = tmp_seaice(i)
822               tmp_tslab(i) = RTT-1.8               tmp_tslab(i) = RTT-1.8
823            ENDIF            ENDIF
824            !  
825            ! faire fondre de la glace si temperature est superieure a 0:            ! faire fondre de la glace si temperature est superieure a 0:
826            !  
827            IF ((tmp_tslab(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN            IF ((tmp_tslab(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN
828               zz = cyang/cbing * (tmp_tslab(i)-RTT)               zz = cyang/cbing * (tmp_tslab(i)-RTT)
829               zz = MIN(zz,tmp_seaice(i))               zz = MIN(zz, tmp_seaice(i))
830               tmp_seaice(i) = tmp_seaice(i) - zz               tmp_seaice(i) = tmp_seaice(i) - zz
831               seaice(i) = tmp_seaice(i)               seaice(i) = tmp_seaice(i)
832               tmp_tslab(i) = tmp_tslab(i) - zz*cbing/cyang               tmp_tslab(i) = tmp_tslab(i) - zz*cbing/cyang
833            ENDIF            ENDIF
834            !  
835            ! limiter la glace de mer a 10 metres (10000 kg/m2)            ! limiter la glace de mer a 10 metres (10000 kg/m2)
836            !  
837            IF(tmp_seaice(i).GT.45.) THEN            IF(tmp_seaice(i).GT.45.) THEN
838               tmp_seaice(i) = MIN(tmp_seaice(i),10000.0)               tmp_seaice(i) = MIN(tmp_seaice(i), 10000.0)
839            ELSE            ELSE
840               tmp_seaice(i) = 0.               tmp_seaice(i) = 0.
841            ENDIF            ENDIF
842            seaice(i) = tmp_seaice(i)            seaice(i) = tmp_seaice(i)
843            siceh(i)=tmp_seaice(i)/1000. !en metres            siceh(i)=tmp_seaice(i)/1000. !en metres
844            !  
845            ! determiner la nature du sol (glace de mer ou ocean libre):            ! determiner la nature du sol (glace de mer ou ocean libre):
846            !  
847            ! on fait dependre la fraction de seaice "pctsrf(i,is_sic)"            ! on fait dependre la fraction de seaice "pctsrf(i, is_sic)"
848            ! de l'epaisseur de seaice :            ! de l'epaisseur de seaice :
849            ! pctsrf(i,is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm            ! pctsrf(i, is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm
850            ! et pctsrf(i,is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur            ! et pctsrf(i, is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur
851            !  
852            pctsrf_slab(i,is_sic)=MIN(siceh(i)/0.20, &            pctsrf_slab(i, is_sic)=MIN(siceh(i)/0.20, &
853                 &                      1.-(pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)))                 1.-(pctsrf_slab(i, is_ter)+pctsrf_slab(i, is_lic)))
854            pctsrf_slab(i,is_oce)=1.0 - &            pctsrf_slab(i, is_oce)=1.0 - &
855                 &      (pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)+pctsrf_slab(i,is_sic))                 (pctsrf_slab(i, is_ter)+pctsrf_slab(i, is_lic)+pctsrf_slab(i, is_sic))
856         ENDIF !pctsrf         ENDIF !pctsrf
857      ENDDO      ENDDO
858      !  
859      ! Calculer le bilan du flux de chaleur au sol :      ! Calculer le bilan du flux de chaleur au sol :
860      !  
861      DO i = 1, klon      DO i = 1, klon
862         za = radsol(i) + fluxo(i)         za = radsol(i) + fluxo(i)
863         zb = fluxg(i)         zb = fluxg(i)
864         IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &         IF((pctsrf_slab(i, is_oce).GT.epsfra).OR. &
865              &   (pctsrf_slab(i,is_sic).GT.epsfra)) THEN              (pctsrf_slab(i, is_sic).GT.epsfra)) THEN
866            slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i,is_oce) &            slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i, is_oce) &
867                 &             +zb*pctsrf_slab(i,is_sic))/ FLOAT(lmt_pas)                 +zb*pctsrf_slab(i, is_sic))/ FLOAT(lmt_pas)
868         ENDIF         ENDIF
869      ENDDO !klon      ENDDO !klon
870      !  
871      ! calcul tslab      ! calcul tslab
872      !  
873      IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee      IF (MOD(itap, lmt_pas).EQ.0) THEN !fin de journee
874         DO i = 1, klon         DO i = 1, klon
875            IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &            IF ((pctsrf_slab(i, is_oce).GT.epsfra).OR. &
876                 &    (pctsrf_slab(i,is_sic).GT.epsfra)) THEN                 (pctsrf_slab(i, is_sic).GT.epsfra)) THEN
877               tmp_tslab(i) = tmp_tslab(i) + &               tmp_tslab(i) = tmp_tslab(i) + &
878                    & (slab_bils(i)-lmt_bils(i)) &                    (slab_bils(i)-lmt_bils(i)) &
879                    &                         /cyang*unjour                    /cyang*unjour
880               ! on remet l'accumulation a 0               ! on remet l'accumulation a 0
881               slab_bils(i) = 0.               slab_bils(i) = 0.
882            ENDIF !pctsrf            ENDIF !pctsrf
883         ENDDO !klon         ENDDO !klon
884      ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN      ENDIF !(MOD(itap, lmt_pas).EQ.0) THEN
885      !  
886      tslab = tmp_tslab      tslab = tmp_tslab
887      seaice = tmp_seaice      seaice = tmp_seaice
888    END SUBROUTINE interfoce_slab    END SUBROUTINE interfoce_slab
889    
890    !************************    !************************
891    
892    SUBROUTINE interfoce_lim(itime, dtime, jour, &    SUBROUTINE interfoce_lim(itime, dtime, jour,  &
893         & klon, nisurf, knon, knindex, &         klon, nisurf, knon, knindex,  &
894         & debut,  &         debut,  &
895         & lmt_sst, pctsrf_new)         lmt_sst, pctsrf_new)
896    
897      ! Cette routine sert d'interface entre le modele atmospherique et un fichier      ! Cette routine sert d'interface entre le modele atmospherique et
898      ! de conditions aux limites      ! un fichier de conditions aux limites
899      !  
900      ! L. Fairhead 02/2000      ! L. Fairhead 02/2000
     !  
     ! input:  
     !   itime        numero du pas de temps courant  
     !   dtime        pas de temps de la physique (en s)  
     !   jour         jour a lire dans l'annee  
     !   nisurf       index de la surface a traiter (1 = sol continental)  
     !   knon         nombre de points dans le domaine a traiter  
     !   knindex      index des points de la surface a traiter  
     !   klon         taille de la grille  
     !   debut        logical: 1er appel a la physique (initialisation)  
     !  
     ! output:  
     !   lmt_sst      SST lues dans le fichier de CL  
     !   pctsrf_new   sous-maille fractionnelle  
     !  
901    
902      use abort_gcm_m, only: abort_gcm      use abort_gcm_m, only: abort_gcm
903      use indicesol      use indicesol
904    
905      ! Parametres d'entree      integer, intent(IN) :: itime ! numero du pas de temps courant
906      integer, intent(IN) :: itime      real , intent(IN) :: dtime ! pas de temps de la physique (en s)
907      real   , intent(IN) :: dtime      integer, intent(IN) :: jour ! jour a lire dans l'annee
908      integer, intent(IN) :: jour      integer, intent(IN) :: nisurf ! index de la surface a traiter (1 = sol continental)
909      integer, intent(IN) :: nisurf      integer, intent(IN) :: knon ! nombre de points dans le domaine a traiter
910      integer, intent(IN) :: knon      integer, intent(IN) :: klon ! taille de la grille
911      integer, intent(IN) :: klon      integer, dimension(klon), intent(in) :: knindex ! index des points de la surface a traiter
912      integer, dimension(klon), intent(in) :: knindex      logical, intent(IN) :: debut ! logical: 1er appel a la physique (initialisation)
     logical, intent(IN) :: debut  
913    
914      ! Parametres de sortie      ! Parametres de sortie
915        ! output:
916        ! lmt_sst SST lues dans le fichier de CL
917        ! pctsrf_new sous-maille fractionnelle
918      real, intent(out), dimension(klon) :: lmt_sst      real, intent(out), dimension(klon) :: lmt_sst
919      real, intent(out), dimension(klon,nbsrf) :: pctsrf_new      real, intent(out), dimension(klon, nbsrf) :: pctsrf_new
920    
921      ! Variables locales      ! Variables locales
922      integer     :: ii      integer :: ii
923      INTEGER,save :: lmt_pas     ! frequence de lecture des conditions limites      INTEGER, save :: lmt_pas ! frequence de lecture des conditions limites
924      ! (en pas de physique)      ! (en pas de physique)
925      logical,save :: deja_lu    ! pour indiquer que le jour a lire a deja      logical, save :: deja_lu ! pour indiquer que le jour a lire a deja
926      ! lu pour une surface precedente      ! lu pour une surface precedente
927      integer,save :: jour_lu      integer, save :: jour_lu
928      integer      :: ierr      integer :: ierr
929      character (len = 20) :: modname = 'interfoce_lim'      character (len = 20) :: modname = 'interfoce_lim'
930      character (len = 80) :: abort_message      character (len = 80) :: abort_message
931      logical, save     :: newlmt = .TRUE.      logical, save :: newlmt = .TRUE.
932      logical, save     :: check = .FALSE.      logical, save :: check = .FALSE.
933      ! Champs lus dans le fichier de CL      ! Champs lus dans le fichier de CL
934      real, allocatable , save, dimension(:) :: sst_lu, rug_lu, nat_lu      real, allocatable , save, dimension(:) :: sst_lu, rug_lu, nat_lu
935      real, allocatable , save, dimension(:,:) :: pct_tmp      real, allocatable , save, dimension(:, :) :: pct_tmp
936      !  
937      ! quelques variables pour netcdf      ! quelques variables pour netcdf
938      !  
939      include "netcdf.inc"      include "netcdf.inc"
940      integer              :: nid, nvarid      integer :: nid, nvarid
941      integer, dimension(2) :: start, epais      integer, dimension(2) :: start, epais
942      !  
943      ! Fin déclaration      ! --------------------------------------------------
     !  
944    
945      if (debut .and. .not. allocated(sst_lu)) then      if (debut .and. .not. allocated(sst_lu)) then
946         lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour         lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
947         jour_lu = jour - 1         jour_lu = jour - 1
948         allocate(sst_lu(klon))         allocate(sst_lu(klon))
949         allocate(nat_lu(klon))         allocate(nat_lu(klon))
950         allocate(pct_tmp(klon,nbsrf))         allocate(pct_tmp(klon, nbsrf))
951      endif      endif
952    
953      if ((jour - jour_lu) /= 0) deja_lu = .false.      if ((jour - jour_lu) /= 0) deja_lu = .false.
954    
955      if (check) write(*,*)modname,' :: jour, jour_lu, deja_lu', jour, jour_lu, &      if (check) write(*, *)modname, ' :: jour, jour_lu, deja_lu', jour, jour_lu, &
956           deja_lu           deja_lu
957      if (check) write(*,*)modname,' :: itime, lmt_pas ', itime, lmt_pas,dtime      if (check) write(*, *)modname, ' :: itime, lmt_pas ', itime, lmt_pas, dtime
958    
959      ! Tester d'abord si c'est le moment de lire le fichier      ! Tester d'abord si c'est le moment de lire le fichier
960      if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then      if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then
961         !  
962         ! Ouverture du fichier         ! Ouverture du fichier
963         !  
964         ierr = NF_OPEN ('limit.nc', NF_NOWRITE,nid)         ierr = NF_OPEN ('limit.nc', NF_NOWRITE, nid)
965         if (ierr.NE.NF_NOERR) then         if (ierr.NE.NF_NOERR) then
966            abort_message &            abort_message &
967                 = 'Pb d''ouverture du fichier de conditions aux limites'                 = 'Pb d''ouverture du fichier de conditions aux limites'
968            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
969         endif         endif
970         !  
971         ! La tranche de donnees a lire:         ! La tranche de donnees a lire:
972         !  
973         start(1) = 1         start(1) = 1
974         start(2) = jour         start(2) = jour
975         epais(1) = klon         epais(1) = klon
976         epais(2) = 1         epais(2) = 1
977         !  
978         if (newlmt) then         if (newlmt) then
979            !  
980            ! Fraction "ocean"            ! Fraction "ocean"
981            !  
982            ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)            ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)
983            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
984               abort_message = 'Le champ <FOCE> est absent'               abort_message = 'Le champ <FOCE> est absent'
985               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
986            endif            endif
987            ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_oce))            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_oce))
988            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
989               abort_message = 'Lecture echouee pour <FOCE>'               abort_message = 'Lecture echouee pour <FOCE>'
990               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
991            endif            endif
992            !  
993            ! Fraction "glace de mer"            ! Fraction "glace de mer"
994            !  
995            ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)            ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)
996            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
997               abort_message = 'Le champ <FSIC> est absent'               abort_message = 'Le champ <FSIC> est absent'
998               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
999            endif            endif
1000            ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_sic))            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_sic))
1001            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1002               abort_message = 'Lecture echouee pour <FSIC>'               abort_message = 'Lecture echouee pour <FSIC>'
1003               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1004            endif            endif
1005            !  
1006            ! Fraction "terre"            ! Fraction "terre"
1007            !  
1008            ierr = NF_INQ_VARID(nid, 'FTER', nvarid)            ierr = NF_INQ_VARID(nid, 'FTER', nvarid)
1009            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1010               abort_message = 'Le champ <FTER> est absent'               abort_message = 'Le champ <FTER> est absent'
1011               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1012            endif            endif
1013            ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_ter))            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_ter))
1014            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1015               abort_message = 'Lecture echouee pour <FTER>'               abort_message = 'Lecture echouee pour <FTER>'
1016               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1017            endif            endif
1018            !  
1019            ! Fraction "glacier terre"            ! Fraction "glacier terre"
1020            !  
1021            ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)            ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)
1022            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1023               abort_message = 'Le champ <FLIC> est absent'               abort_message = 'Le champ <FLIC> est absent'
1024               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1025            endif            endif
1026            ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais,pct_tmp(1,is_lic))            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_lic))
1027            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1028               abort_message = 'Lecture echouee pour <FLIC>'               abort_message = 'Lecture echouee pour <FLIC>'
1029               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1030            endif            endif
1031            !  
1032         else  ! on en est toujours a rnatur         else ! on en est toujours a rnatur
1033            !  
1034            ierr = NF_INQ_VARID(nid, 'NAT', nvarid)            ierr = NF_INQ_VARID(nid, 'NAT', nvarid)
1035            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1036               abort_message = 'Le champ <NAT> est absent'               abort_message = 'Le champ <NAT> est absent'
1037               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1038            endif            endif
1039            ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, nat_lu)            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, nat_lu)
1040            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1041               abort_message = 'Lecture echouee pour <NAT>'               abort_message = 'Lecture echouee pour <NAT>'
1042               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1043            endif            endif
1044            !  
1045            ! Remplissage des fractions de surface            ! Remplissage des fractions de surface
1046            ! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice            ! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice
1047            !  
1048            pct_tmp = 0.0            pct_tmp = 0.0
1049            do ii = 1, klon            do ii = 1, klon
1050               pct_tmp(ii,nint(nat_lu(ii)) + 1) = 1.               pct_tmp(ii, nint(nat_lu(ii)) + 1) = 1.
1051            enddo            enddo
1052    
1053            !  
1054            !  On se retrouve avec ocean en 1 et terre en 2 alors qu'on veut le contraire            ! On se retrouve avec ocean en 1 et terre en 2 alors qu'on veut le contraire
1055            !  
1056            pctsrf_new = pct_tmp            pctsrf_new = pct_tmp
1057            pctsrf_new (:,2)= pct_tmp (:,1)            pctsrf_new (:, 2)= pct_tmp (:, 1)
1058            pctsrf_new (:,1)= pct_tmp (:,2)            pctsrf_new (:, 1)= pct_tmp (:, 2)
1059            pct_tmp = pctsrf_new            pct_tmp = pctsrf_new
1060         endif ! fin test sur newlmt         endif ! fin test sur newlmt
1061         !  
1062         ! Lecture SST         ! Lecture SST
1063         !  
1064         ierr = NF_INQ_VARID(nid, 'SST', nvarid)         ierr = NF_INQ_VARID(nid, 'SST', nvarid)
1065         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1066            abort_message = 'Le champ <SST> est absent'            abort_message = 'Le champ <SST> est absent'
1067            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1068         endif         endif
1069         ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, sst_lu)         ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, sst_lu)
1070         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1071            abort_message = 'Lecture echouee pour <SST>'            abort_message = 'Lecture echouee pour <SST>'
1072            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1073         endif         endif
1074    
1075         !  
1076         ! Fin de lecture         ! Fin de lecture
1077         !  
1078         ierr = NF_CLOSE(nid)         ierr = NF_CLOSE(nid)
1079         deja_lu = .true.         deja_lu = .true.
1080         jour_lu = jour         jour_lu = jour
1081      endif      endif
1082      !  
1083      ! Recopie des variables dans les champs de sortie      ! Recopie des variables dans les champs de sortie
1084      !  
1085      lmt_sst = 999999999.      lmt_sst = 999999999.
1086      do ii = 1, knon      do ii = 1, knon
1087         lmt_sst(ii) = sst_lu(knindex(ii))         lmt_sst(ii) = sst_lu(knindex(ii))
1088      enddo      enddo
1089    
1090      pctsrf_new(:,is_oce) = pct_tmp(:,is_oce)      pctsrf_new(:, is_oce) = pct_tmp(:, is_oce)
1091      pctsrf_new(:,is_sic) = pct_tmp(:,is_sic)      pctsrf_new(:, is_sic) = pct_tmp(:, is_sic)
1092    
1093    END SUBROUTINE interfoce_lim    END SUBROUTINE interfoce_lim
1094    
1095    !************************    !************************
1096    
1097    SUBROUTINE interfsur_lim(itime, dtime, jour, &    SUBROUTINE interfsur_lim(itime, dtime, jour,  &
1098         & klon, nisurf, knon, knindex, &         klon, nisurf, knon, knindex,  &
1099         & debut,  &         debut,  &
1100         & lmt_alb, lmt_rug)         lmt_alb, lmt_rug)
1101    
1102      ! Cette routine sert d'interface entre le modèle atmosphérique et      ! Cette routine sert d'interface entre le modèle atmosphérique et
1103      ! un fichier de conditions aux limites.      ! un fichier de conditions aux limites.
1104      !  
1105      ! L. Fairhead 02/2000      ! L. Fairhead 02/2000
     !  
     ! input:  
     !   itime        numero du pas de temps courant  
     !   dtime        pas de temps de la physique (en s)  
     !   jour         jour a lire dans l'annee  
     !   nisurf       index de la surface a traiter (1 = sol continental)  
     !   knon         nombre de points dans le domaine a traiter  
     !   knindex      index des points de la surface a traiter  
     !   klon         taille de la grille  
     !   debut        logical: 1er appel a la physique (initialisation)  
     !  
     ! output:  
     !   lmt_sst      SST lues dans le fichier de CL  
     !   lmt_alb      Albedo lu  
     !   lmt_rug      longueur de rugosité lue  
     !   pctsrf_new   sous-maille fractionnelle  
     !  
1106    
1107      use abort_gcm_m, only: abort_gcm      use abort_gcm_m, only: abort_gcm
1108    
1109      ! Parametres d'entree      ! Parametres d'entree
1110        ! input:
1111        ! itime numero du pas de temps courant
1112        ! dtime pas de temps de la physique (en s)
1113        ! jour jour a lire dans l'annee
1114        ! nisurf index de la surface a traiter (1 = sol continental)
1115        ! knon nombre de points dans le domaine a traiter
1116        ! knindex index des points de la surface a traiter
1117        ! klon taille de la grille
1118        ! debut logical: 1er appel a la physique (initialisation)
1119      integer, intent(IN) :: itime      integer, intent(IN) :: itime
1120      real   , intent(IN) :: dtime      real , intent(IN) :: dtime
1121      integer, intent(IN) :: jour      integer, intent(IN) :: jour
1122      integer, intent(IN) :: nisurf      integer, intent(IN) :: nisurf
1123      integer, intent(IN) :: knon      integer, intent(IN) :: knon
# Line 1869  CONTAINS Line 1126  CONTAINS
1126      logical, intent(IN) :: debut      logical, intent(IN) :: debut
1127    
1128      ! Parametres de sortie      ! Parametres de sortie
1129        ! output:
1130        ! lmt_sst SST lues dans le fichier de CL
1131        ! lmt_alb Albedo lu
1132        ! lmt_rug longueur de rugosité lue
1133        ! pctsrf_new sous-maille fractionnelle
1134      real, intent(out), dimension(klon) :: lmt_alb      real, intent(out), dimension(klon) :: lmt_alb
1135      real, intent(out), dimension(klon) :: lmt_rug      real, intent(out), dimension(klon) :: lmt_rug
1136    
1137      ! Variables locales      ! Variables locales
1138      integer     :: ii      integer :: ii
1139      integer,save :: lmt_pas     ! frequence de lecture des conditions limites      integer, save :: lmt_pas ! frequence de lecture des conditions limites
1140      ! (en pas de physique)      ! (en pas de physique)
1141      logical,save :: deja_lu_sur! pour indiquer que le jour a lire a deja      logical, save :: deja_lu_sur! pour indiquer que le jour a lire a deja
1142      ! lu pour une surface precedente      ! lu pour une surface precedente
1143      integer,save :: jour_lu_sur      integer, save :: jour_lu_sur
1144      integer      :: ierr      integer :: ierr
1145      character (len = 20) :: modname = 'interfsur_lim'      character (len = 20) :: modname = 'interfsur_lim'
1146      character (len = 80) :: abort_message      character (len = 80) :: abort_message
1147      logical,save     :: newlmt = .false.      logical, save :: newlmt = .false.
1148      logical,save     :: check = .false.      logical, save :: check = .false.
1149      ! Champs lus dans le fichier de CL      ! Champs lus dans le fichier de CL
1150      real, allocatable , save, dimension(:) :: alb_lu, rug_lu      real, allocatable , save, dimension(:) :: alb_lu, rug_lu
1151      !  
1152      ! quelques variables pour netcdf      ! quelques variables pour netcdf
1153      !  
1154      include "netcdf.inc"      include "netcdf.inc"
1155      integer ,save             :: nid, nvarid      integer , save :: nid, nvarid
1156      integer, dimension(2),save :: start, epais      integer, dimension(2), save :: start, epais
1157      !  
1158      ! Fin déclaration      !------------------------------------------------------------
     !  
1159    
1160      if (debut) then      if (debut) then
1161         lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour         lmt_pas = nint(86400./dtime * 1.0) ! pour une lecture une fois par jour
# Line 1905  CONTAINS Line 1166  CONTAINS
1166    
1167      if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false.      if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false.
1168    
1169      if (check) write(*,*)modname,':: jour_lu_sur, deja_lu_sur', jour_lu_sur, &      if (check) write(*, *)modname, ':: jour_lu_sur, deja_lu_sur', jour_lu_sur, &
1170           deja_lu_sur           deja_lu_sur
1171      if (check) write(*,*)modname,':: itime, lmt_pas', itime, lmt_pas      if (check) write(*, *)modname, ':: itime, lmt_pas', itime, lmt_pas
1172    
1173      ! Tester d'abord si c'est le moment de lire le fichier      ! Tester d'abord si c'est le moment de lire le fichier
1174      if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu_sur) then      if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu_sur) then
1175         !  
1176         ! Ouverture du fichier         ! Ouverture du fichier
1177         !  
1178         ierr = NF_OPEN ('limit.nc', NF_NOWRITE,nid)         ierr = NF_OPEN ('limit.nc', NF_NOWRITE, nid)
1179         if (ierr.NE.NF_NOERR) then         if (ierr.NE.NF_NOERR) then
1180            abort_message &            abort_message &
1181                 = 'Pb d''ouverture du fichier de conditions aux limites'                 = 'Pb d''ouverture du fichier de conditions aux limites'
1182            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1183         endif         endif
1184         !  
1185         ! La tranche de donnees a lire:         ! La tranche de donnees a lire:
1186    
1187         start(1) = 1         start(1) = 1
1188         start(2) = jour         start(2) = jour
1189         epais(1) = klon         epais(1) = klon
1190         epais(2) = 1         epais(2) = 1
1191         !  
1192         ! Lecture Albedo         ! Lecture Albedo
1193         !  
1194         ierr = NF_INQ_VARID(nid, 'ALB', nvarid)         ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
1195         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1196            abort_message = 'Le champ <ALB> est absent'            abort_message = 'Le champ <ALB> est absent'
1197            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1198         endif         endif
1199         ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, alb_lu)         ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, alb_lu)
1200         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1201            abort_message = 'Lecture echouee pour <ALB>'            abort_message = 'Lecture echouee pour <ALB>'
1202            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1203         endif         endif
1204         !  
1205         ! Lecture rugosité         ! Lecture rugosité
1206         !  
1207         ierr = NF_INQ_VARID(nid, 'RUG', nvarid)         ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
1208         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1209            abort_message = 'Le champ <RUG> est absent'            abort_message = 'Le champ <RUG> est absent'
1210            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1211         endif         endif
1212         ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)         ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, rug_lu)
1213         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1214            abort_message = 'Lecture echouee pour <RUG>'            abort_message = 'Lecture echouee pour <RUG>'
1215            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1216         endif         endif
1217    
1218         !  
1219         ! Fin de lecture         ! Fin de lecture
1220         !  
1221         ierr = NF_CLOSE(nid)         ierr = NF_CLOSE(nid)
1222         deja_lu_sur = .true.         deja_lu_sur = .true.
1223         jour_lu_sur = jour         jour_lu_sur = jour
1224      endif      endif
1225      !  
1226      ! Recopie des variables dans les champs de sortie      ! Recopie des variables dans les champs de sortie
1227      !  
1228  !!$  lmt_alb(:) = 0.0  !!$ lmt_alb = 0.0
1229  !!$  lmt_rug(:) = 0.0  !!$ lmt_rug = 0.0
1230      lmt_alb(:) = 999999.      lmt_alb = 999999.
1231      lmt_rug(:) = 999999.      lmt_rug = 999999.
1232      DO ii = 1, knon      DO ii = 1, knon
1233         lmt_alb(ii) = alb_lu(knindex(ii))         lmt_alb(ii) = alb_lu(knindex(ii))
1234         lmt_rug(ii) = rug_lu(knindex(ii))         lmt_rug(ii) = rug_lu(knindex(ii))
# Line 1977  CONTAINS Line 1238  CONTAINS
1238    
1239    !************************    !************************
1240    
1241    SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime, &    SUBROUTINE calcul_fluxs( klon, knon, nisurf, dtime,  &
1242         & tsurf, p1lay, cal, beta, coef1lay, ps, &         tsurf, p1lay, cal, beta, coef1lay, ps,  &
1243         & precip_rain, precip_snow, snow, qsurf, &         precip_rain, precip_snow, snow, qsurf,  &
1244         & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &         radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay,  &
1245         & petAcoef, peqAcoef, petBcoef, peqBcoef, &         petAcoef, peqAcoef, petBcoef, peqBcoef,  &
1246         & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
1247    
1248      ! Cette routine calcule les fluxs en h et q a l'interface et eventuellement      ! Cette routine calcule les fluxs en h et q a l'interface et eventuellement
1249      ! une temperature de surface (au cas ou ok_veget = false)      ! une temperature de surface (au cas ou ok_veget = false)
1250      !  
1251      ! L. Fairhead 4/2000      ! L. Fairhead 4/2000
1252      !  
1253      ! input:      ! input:
1254      !   knon         nombre de points a traiter      ! knon nombre de points a traiter
1255      !   nisurf       surface a traiter      ! nisurf surface a traiter
1256      !   tsurf        temperature de surface      ! tsurf temperature de surface
1257      !   p1lay        pression 1er niveau (milieu de couche)      ! p1lay pression 1er niveau (milieu de couche)
1258      !   cal          capacite calorifique du sol      ! cal capacite calorifique du sol
1259      !   beta         evap reelle      ! beta evap reelle
1260      !   coef1lay     coefficient d'echange      ! coef1lay coefficient d'echange
1261      !   ps           pression au sol      ! ps pression au sol
1262      !   precip_rain  precipitations liquides      ! precip_rain precipitations liquides
1263      !   precip_snow  precipitations solides      ! precip_snow precipitations solides
1264      !   snow         champs hauteur de neige      ! snow champs hauteur de neige
1265      !   runoff       runoff en cas de trop plein      ! runoff runoff en cas de trop plein
1266      !   petAcoef     coeff. A de la resolution de la CL pour t      ! petAcoef coeff. A de la resolution de la CL pour t
1267      !   peqAcoef     coeff. A de la resolution de la CL pour q      ! peqAcoef coeff. A de la resolution de la CL pour q
1268      !   petBcoef     coeff. B de la resolution de la CL pour t      ! petBcoef coeff. B de la resolution de la CL pour t
1269      !   peqBcoef     coeff. B de la resolution de la CL pour q      ! peqBcoef coeff. B de la resolution de la CL pour q
1270      !   radsol       rayonnement net aus sol (LW + SW)      ! radsol rayonnement net aus sol (LW + SW)
1271      !   dif_grnd     coeff. diffusion vers le sol profond      ! dif_grnd coeff. diffusion vers le sol profond
1272      !  
1273      ! output:      ! output:
1274      !   tsurf_new    temperature au sol      ! tsurf_new temperature au sol
1275      !   qsurf        humidite de l'air au dessus du sol      ! qsurf humidite de l'air au dessus du sol
1276      !   fluxsens     flux de chaleur sensible      ! fluxsens flux de chaleur sensible
1277      !   fluxlat      flux de chaleur latente      ! fluxlat flux de chaleur latente
1278      !   dflux_s      derivee du flux de chaleur sensible / Ts      ! dflux_s derivee du flux de chaleur sensible / Ts
1279      !   dflux_l      derivee du flux de chaleur latente  / Ts      ! dflux_l derivee du flux de chaleur latente / Ts
1280      !  
1281    
1282      use indicesol      use indicesol
1283      use abort_gcm_m, only: abort_gcm      use abort_gcm_m, only: abort_gcm
1284      use yoethf      use yoethf_m
1285      use fcttre, only: thermcep, foeew, qsats, qsatl, foede, dqsats, dqsatl      use fcttre, only: thermcep, foeew, qsats, qsatl, foede, dqsats, dqsatl
1286      use YOMCST      use SUPHEC_M
1287    
1288      ! Parametres d'entree      ! Parametres d'entree
1289      integer, intent(IN) :: knon, nisurf, klon      integer, intent(IN) :: knon, nisurf, klon
1290      real   , intent(IN) :: dtime      real , intent(IN) :: dtime
1291      real, dimension(klon), intent(IN) :: petAcoef, peqAcoef      real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
1292      real, dimension(klon), intent(IN) :: petBcoef, peqBcoef      real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
1293      real, dimension(klon), intent(IN) :: ps, q1lay      real, dimension(klon), intent(IN) :: ps, q1lay
# Line 2047  CONTAINS Line 1308  CONTAINS
1308      real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef      real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
1309      real, dimension(klon) :: zx_sl, zx_k1      real, dimension(klon) :: zx_sl, zx_k1
1310      real, dimension(klon) :: zx_q_0 , d_ts      real, dimension(klon) :: zx_q_0 , d_ts
1311      real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh      real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
1312      real                  :: bilan_f, fq_fonte      real :: bilan_f, fq_fonte
1313      REAL                  :: subli, fsno      REAL :: subli, fsno
1314      REAL                  :: qsat_new, q1_new      REAL :: qsat_new, q1_new
1315      real, parameter :: t_grnd = 271.35, t_coup = 273.15      real, parameter :: t_grnd = 271.35, t_coup = 273.15
1316      !! PB temporaire en attendant mieux pour le modele de neige      !! PB temporaire en attendant mieux pour le modele de neige
1317      REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)      REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
1318      !  
1319      logical, save         :: check = .false.      logical, save :: check = .false.
1320      character (len = 20)  :: modname = 'calcul_fluxs'      character (len = 20) :: modname = 'calcul_fluxs'
1321      logical, save         :: fonte_neige = .false.      logical, save :: fonte_neige = .false.
1322      real, save            :: max_eau_sol = 150.0      real, save :: max_eau_sol = 150.0
1323      character (len = 80) :: abort_message      character (len = 80) :: abort_message
1324      logical,save         :: first = .true.,second=.false.      logical, save :: first = .true., second=.false.
1325    
1326      if (check) write(*,*)'Entree ', modname,' surface = ',nisurf      if (check) write(*, *)'Entree ', modname, ' surface = ', nisurf
1327    
1328      IF (check) THEN      IF (check) THEN
1329         WRITE(*,*)' radsol (min, max)' &         WRITE(*, *)' radsol (min, max)' &
1330              &     , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))              , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
1331         !!CALL flush(6)         !!CALL flush(6)
1332      ENDIF      ENDIF
1333    
1334      if (size(coastalflow) /= knon .AND. nisurf == is_ter) then      if (size(coastalflow) /= knon .AND. nisurf == is_ter) then
1335         write(*,*)'Bizarre, le nombre de points continentaux'         write(*, *)'Bizarre, le nombre de points continentaux'
1336         write(*,*)'a change entre deux appels. J''arrete ...'         write(*, *)'a change entre deux appels. J''arrete ...'
1337         abort_message='Pb run_off'         abort_message='Pb run_off'
1338         call abort_gcm(modname,abort_message,1)         call abort_gcm(modname, abort_message, 1)
1339      endif      endif
1340      !  
1341      ! Traitement neige et humidite du sol      ! Traitement neige et humidite du sol
     !  
 !!$  WRITE(*,*)'test calcul_flux, surface ', nisurf  
     !!PB test  
 !!$    if (nisurf == is_oce) then  
 !!$      snow = 0.  
 !!$      qsol = max_eau_sol  
 !!$    else  
 !!$      where (precip_snow > 0.) snow = snow + (precip_snow * dtime)  
 !!$      where (snow > epsilon(snow)) snow = max(0.0, snow - (evap * dtime))  
 !!$!      snow = max(0.0, snow + (precip_snow - evap) * dtime)  
 !!$      where (precip_rain > 0.) qsol = qsol + (precip_rain - evap) * dtime  
 !!$    endif  
 !!$    IF (nisurf /= is_ter) qsol = max_eau_sol  
1342    
     !  
1343      ! Initialisation      ! Initialisation
1344      !  
1345      evap = 0.      evap = 0.
1346      fluxsens=0.      fluxsens=0.
1347      fluxlat=0.      fluxlat=0.
1348      dflux_s = 0.      dflux_s = 0.
1349      dflux_l = 0.      dflux_l = 0.
1350      !  
1351      ! zx_qs = qsat en kg/kg      ! zx_qs = qsat en kg/kg
1352      !  
1353      DO i = 1, knon      DO i = 1, knon
1354         zx_pkh(i) = (ps(i)/ps(i))**RKAPPA         zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
1355         IF (thermcep) THEN         IF (thermcep) THEN
1356            zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))            zdelta=MAX(0., SIGN(1., rtt-tsurf(i)))
1357            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
1358            zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))            zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
1359            zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)            zx_qs= r2es * FOEEW(tsurf(i), zdelta)/ps(i)
1360            zx_qs=MIN(0.5,zx_qs)            zx_qs=MIN(0.5, zx_qs)
1361            zcor=1./(1.-retv*zx_qs)            zcor=1./(1.-retv*zx_qs)
1362            zx_qs=zx_qs*zcor            zx_qs=zx_qs*zcor
1363            zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &            zx_dq_s_dh = FOEDE(tsurf(i), zdelta, zcvm5, zx_qs, zcor) &
1364                 &                 /RLVTT / zx_pkh(i)                 /RLVTT / zx_pkh(i)
1365         ELSE         ELSE
1366            IF (tsurf(i).LT.t_coup) THEN            IF (tsurf(i).LT.t_coup) THEN
1367               zx_qs = qsats(tsurf(i)) / ps(i)               zx_qs = qsats(tsurf(i)) / ps(i)
1368               zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &               zx_dq_s_dh = dqsats(tsurf(i), zx_qs)/RLVTT &
1369                    &                    / zx_pkh(i)                    / zx_pkh(i)
1370            ELSE            ELSE
1371               zx_qs = qsatl(tsurf(i)) / ps(i)               zx_qs = qsatl(tsurf(i)) / ps(i)
1372               zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &               zx_dq_s_dh = dqsatl(tsurf(i), zx_qs)/RLVTT &
1373                    &               / zx_pkh(i)                    / zx_pkh(i)
1374            ENDIF            ENDIF
1375         ENDIF         ENDIF
1376         zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh         zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
1377         zx_qsat(i) = zx_qs         zx_qsat(i) = zx_qs
1378         zx_coef(i) = coef1lay(i) &         zx_coef(i) = coef1lay(i) &
1379              & * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &              * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
1380              & * p1lay(i)/(RD*t1lay(i))              * p1lay(i)/(RD*t1lay(i))
1381    
1382      ENDDO      ENDDO
1383    
1384      ! === Calcul de la temperature de surface ===      ! === Calcul de la temperature de surface ===
1385      !  
1386      ! zx_sl = chaleur latente d'evaporation ou de sublimation      ! zx_sl = chaleur latente d'evaporation ou de sublimation
1387      !  
1388      do i = 1, knon      do i = 1, knon
1389         zx_sl(i) = RLVTT         zx_sl(i) = RLVTT
1390         if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT         if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
# Line 2148  CONTAINS Line 1395  CONTAINS
1395         ! Q         ! Q
1396         zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)         zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
1397         zx_mq(i) = beta(i) * zx_k1(i) * &         zx_mq(i) = beta(i) * zx_k1(i) * &
1398              &             (peqAcoef(i) - zx_qsat(i) &              (peqAcoef(i) - zx_qsat(i) &
1399              &                          + zx_dq_s_dt(i) * tsurf(i)) &              + zx_dq_s_dt(i) * tsurf(i)) &
1400              &             / zx_oq(i)              / zx_oq(i)
1401         zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &         zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
1402              &                              / zx_oq(i)              / zx_oq(i)
1403    
1404         ! H         ! H
1405         zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)         zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
# Line 2161  CONTAINS Line 1408  CONTAINS
1408    
1409         ! Tsurface         ! Tsurface
1410         tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &         tsurf_new(i) = (tsurf(i) + cal(i)/(RCPD * zx_pkh(i)) * dtime * &
1411              &             (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &              (radsol(i) + zx_mh(i) + zx_sl(i) * zx_mq(i)) &
1412              &                 + dif_grnd(i) * t_grnd * dtime)/ &              + dif_grnd(i) * t_grnd * dtime)/ &
1413              &          ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &              ( 1. - dtime * cal(i)/(RCPD * zx_pkh(i)) * ( &
1414              &                       zx_nh(i) + zx_sl(i) * zx_nq(i)) &                zx_nh(i) + zx_sl(i) * zx_nq(i)) &
1415              &                     + dtime * dif_grnd(i))              + dtime * dif_grnd(i))
1416    
1417    
        !  
1418         ! Y'a-t-il fonte de neige?         ! Y'a-t-il fonte de neige?
1419         !  
1420         !    fonte_neige = (nisurf /= is_oce) .AND. &         ! fonte_neige = (nisurf /= is_oce) .AND. &
1421         !     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &         ! & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
1422         !     & .AND. (tsurf_new(i) >= RTT)         ! & .AND. (tsurf_new(i) >= RTT)
1423         !    if (fonte_neige) tsurf_new(i) = RTT           ! if (fonte_neige) tsurf_new(i) = RTT
1424         d_ts(i) = tsurf_new(i) - tsurf(i)         d_ts(i) = tsurf_new(i) - tsurf(i)
1425         !    zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)         ! zx_h_ts(i) = tsurf_new(i) * RCPD * zx_pkh(i)
1426         !    zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)         ! zx_q_0(i) = zx_qsat(i) + zx_dq_s_dt(i) * d_ts(i)
1427         !== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas         !== flux_q est le flux de vapeur d'eau: kg/(m**2 s) positive vers bas
1428         !== flux_t est le flux de cpt (energie sensible): j/(m**2 s)         !== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
1429         evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)         evap(i) = - zx_mq(i) - zx_nq(i) * tsurf_new(i)
1430         fluxlat(i) = - evap(i) * zx_sl(i)         fluxlat(i) = - evap(i) * zx_sl(i)
# Line 2195  CONTAINS Line 1442  CONTAINS
1442    
1443    !************************    !************************
1444    
1445    SUBROUTINE fonte_neige( klon, knon, nisurf, dtime, &    SUBROUTINE fonte_neige( klon, knon, nisurf, dtime,  &
1446         & tsurf, p1lay, cal, beta, coef1lay, ps, &         tsurf, p1lay, cal, beta, coef1lay, ps,  &
1447         & precip_rain, precip_snow, snow, qsol, &         precip_rain, precip_snow, snow, qsol,  &
1448         & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &         radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay,  &
1449         & petAcoef, peqAcoef, petBcoef, peqBcoef, &         petAcoef, peqAcoef, petBcoef, peqBcoef,  &
1450         & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l,  &
1451         & fqcalving,ffonte,run_off_lic_0)         fqcalving, ffonte, run_off_lic_0)
1452    
1453      ! Routine de traitement de la fonte de la neige dans le cas du traitement      ! Routine de traitement de la fonte de la neige dans le cas du traitement
1454      ! de sol simplifié      ! de sol simplifié
1455      !  
1456      ! LF 03/2001      ! LF 03/2001
1457      ! input:      ! input:
1458      !   knon         nombre de points a traiter      ! knon nombre de points a traiter
1459      !   nisurf       surface a traiter      ! nisurf surface a traiter
1460      !   tsurf        temperature de surface      ! tsurf temperature de surface
1461      !   p1lay        pression 1er niveau (milieu de couche)      ! p1lay pression 1er niveau (milieu de couche)
1462      !   cal          capacite calorifique du sol      ! cal capacite calorifique du sol
1463      !   beta         evap reelle      ! beta evap reelle
1464      !   coef1lay     coefficient d'echange      ! coef1lay coefficient d'echange
1465      !   ps           pression au sol      ! ps pression au sol
1466      !   precip_rain  precipitations liquides      ! precip_rain precipitations liquides
1467      !   precip_snow  precipitations solides      ! precip_snow precipitations solides
1468      !   snow         champs hauteur de neige      ! snow champs hauteur de neige
1469      !   qsol         hauteur d'eau contenu dans le sol      ! qsol hauteur d'eau contenu dans le sol
1470      !   runoff       runoff en cas de trop plein      ! runoff runoff en cas de trop plein
1471      !   petAcoef     coeff. A de la resolution de la CL pour t      ! petAcoef coeff. A de la resolution de la CL pour t
1472      !   peqAcoef     coeff. A de la resolution de la CL pour q      ! peqAcoef coeff. A de la resolution de la CL pour q
1473      !   petBcoef     coeff. B de la resolution de la CL pour t      ! petBcoef coeff. B de la resolution de la CL pour t
1474      !   peqBcoef     coeff. B de la resolution de la CL pour q      ! peqBcoef coeff. B de la resolution de la CL pour q
1475      !   radsol       rayonnement net aus sol (LW + SW)      ! radsol rayonnement net aus sol (LW + SW)
1476      !   dif_grnd     coeff. diffusion vers le sol profond      ! dif_grnd coeff. diffusion vers le sol profond
1477      !  
1478      ! output:      ! output:
1479      !   tsurf_new    temperature au sol      ! tsurf_new temperature au sol
1480      !   fluxsens     flux de chaleur sensible      ! fluxsens flux de chaleur sensible
1481      !   fluxlat      flux de chaleur latente      ! fluxlat flux de chaleur latente
1482      !   dflux_s      derivee du flux de chaleur sensible / Ts      ! dflux_s derivee du flux de chaleur sensible / Ts
1483      !   dflux_l      derivee du flux de chaleur latente  / Ts      ! dflux_l derivee du flux de chaleur latente / Ts
1484      ! in/out:      ! in/out:
1485      !   run_off_lic_0 run off glacier du pas de temps précedent      ! run_off_lic_0 run off glacier du pas de temps précedent
1486      !  
1487    
1488      use indicesol      use indicesol
1489      use YOMCST      use SUPHEC_M
1490      use yoethf      use yoethf_m
1491      use fcttre      use fcttre
1492      !IM cf JLD      !IM cf JLD
1493    
1494      ! Parametres d'entree      ! Parametres d'entree
1495      integer, intent(IN) :: knon, nisurf, klon      integer, intent(IN) :: knon, nisurf, klon
1496      real   , intent(IN) :: dtime      real , intent(IN) :: dtime
1497      real, dimension(klon), intent(IN) :: petAcoef, peqAcoef      real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
1498      real, dimension(klon), intent(IN) :: petBcoef, peqBcoef      real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
1499      real, dimension(klon), intent(IN) :: ps, q1lay      real, dimension(klon), intent(IN) :: ps, q1lay
# Line 2268  CONTAINS Line 1515  CONTAINS
1515      ! Variables locales      ! Variables locales
1516      ! Masse maximum de neige (kg/m2). Au dessus de ce seuil, la neige      ! Masse maximum de neige (kg/m2). Au dessus de ce seuil, la neige
1517      ! en exces "s'ecoule" (calving)      ! en exces "s'ecoule" (calving)
1518      !  real, parameter :: snow_max=1.      ! real, parameter :: snow_max=1.
1519      !IM cf JLD/GK      !IM cf JLD/GK
1520      real, parameter :: snow_max=3000.      real, parameter :: snow_max=3000.
1521      integer :: i      integer :: i
# Line 2277  CONTAINS Line 1524  CONTAINS
1524      real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef      real, dimension(klon) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
1525      real, dimension(klon) :: zx_sl, zx_k1      real, dimension(klon) :: zx_sl, zx_k1
1526      real, dimension(klon) :: zx_q_0 , d_ts      real, dimension(klon) :: zx_q_0 , d_ts
1527      real                  :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh      real :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
1528      real                  :: bilan_f, fq_fonte      real :: bilan_f, fq_fonte
1529      REAL                  :: subli, fsno      REAL :: subli, fsno
1530      REAL, DIMENSION(klon) :: bil_eau_s, snow_evap      REAL, DIMENSION(klon) :: bil_eau_s, snow_evap
1531      real, parameter :: t_grnd = 271.35, t_coup = 273.15      real, parameter :: t_grnd = 271.35, t_coup = 273.15
1532      !! PB temporaire en attendant mieux pour le modele de neige      !! PB temporaire en attendant mieux pour le modele de neige
# Line 2288  CONTAINS Line 1535  CONTAINS
1535      !IM cf JLD/ GKtest      !IM cf JLD/ GKtest
1536      REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)      REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)
1537      ! fin GKtest      ! fin GKtest
1538      !  
1539      logical, save         :: check = .FALSE.      logical, save :: check = .FALSE.
1540      character (len = 20)  :: modname = 'fonte_neige'      character (len = 20) :: modname = 'fonte_neige'
1541      logical, save         :: neige_fond = .false.      logical, save :: neige_fond = .false.
1542      real, save            :: max_eau_sol = 150.0      real, save :: max_eau_sol = 150.0
1543      character (len = 80) :: abort_message      character (len = 80) :: abort_message
1544      logical,save         :: first = .true.,second=.false.      logical, save :: first = .true., second=.false.
1545      real                 :: coeff_rel      real :: coeff_rel
1546    
1547      if (check) write(*,*)'Entree ', modname,' surface = ',nisurf      if (check) write(*, *)'Entree ', modname, ' surface = ', nisurf
1548    
1549      ! Initialisations      ! Initialisations
1550      coeff_rel = dtime/(tau_calv * rday)      coeff_rel = dtime/(tau_calv * rday)
1551      bil_eau_s(:) = 0.      bil_eau_s = 0.
1552      DO i = 1, knon      DO i = 1, knon
1553         zx_pkh(i) = (ps(i)/ps(i))**RKAPPA         zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
1554         IF (thermcep) THEN         IF (thermcep) THEN
1555            zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))            zdelta=MAX(0., SIGN(1., rtt-tsurf(i)))
1556            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
1557            zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))            zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
1558            zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)            zx_qs= r2es * FOEEW(tsurf(i), zdelta)/ps(i)
1559            zx_qs=MIN(0.5,zx_qs)            zx_qs=MIN(0.5, zx_qs)
1560            zcor=1./(1.-retv*zx_qs)            zcor=1./(1.-retv*zx_qs)
1561            zx_qs=zx_qs*zcor            zx_qs=zx_qs*zcor
1562            zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &            zx_dq_s_dh = FOEDE(tsurf(i), zdelta, zcvm5, zx_qs, zcor) &
1563                 &                 /RLVTT / zx_pkh(i)                 /RLVTT / zx_pkh(i)
1564         ELSE         ELSE
1565            IF (tsurf(i).LT.t_coup) THEN            IF (tsurf(i).LT.t_coup) THEN
1566               zx_qs = qsats(tsurf(i)) / ps(i)               zx_qs = qsats(tsurf(i)) / ps(i)
1567               zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &               zx_dq_s_dh = dqsats(tsurf(i), zx_qs)/RLVTT &
1568                    &                    / zx_pkh(i)                    / zx_pkh(i)
1569            ELSE            ELSE
1570               zx_qs = qsatl(tsurf(i)) / ps(i)               zx_qs = qsatl(tsurf(i)) / ps(i)
1571               zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &               zx_dq_s_dh = dqsatl(tsurf(i), zx_qs)/RLVTT &
1572                    &               / zx_pkh(i)                    / zx_pkh(i)
1573            ENDIF            ENDIF
1574         ENDIF         ENDIF
1575         zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh         zx_dq_s_dt(i) = RCPD * zx_pkh(i) * zx_dq_s_dh
1576         zx_qsat(i) = zx_qs         zx_qsat(i) = zx_qs
1577         zx_coef(i) = coef1lay(i) &         zx_coef(i) = coef1lay(i) &
1578              & * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &              * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) &
1579              & * p1lay(i)/(RD*t1lay(i))              * p1lay(i)/(RD*t1lay(i))
1580      ENDDO      ENDDO
1581    
1582      ! === Calcul de la temperature de surface ===      ! === Calcul de la temperature de surface ===
1583      !  
1584      ! zx_sl = chaleur latente d'evaporation ou de sublimation      ! zx_sl = chaleur latente d'evaporation ou de sublimation
1585      !  
1586      do i = 1, knon      do i = 1, knon
1587         zx_sl(i) = RLVTT         zx_sl(i) = RLVTT
1588         if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT         if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
# Line 2346  CONTAINS Line 1593  CONTAINS
1593         ! Q         ! Q
1594         zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)         zx_oq(i) = 1. - (beta(i) * zx_k1(i) * peqBcoef(i) * dtime)
1595         zx_mq(i) = beta(i) * zx_k1(i) * &         zx_mq(i) = beta(i) * zx_k1(i) * &
1596              &             (peqAcoef(i) - zx_qsat(i) &              (peqAcoef(i) - zx_qsat(i) &
1597              &                          + zx_dq_s_dt(i) * tsurf(i)) &              + zx_dq_s_dt(i) * tsurf(i)) &
1598              &             / zx_oq(i)              / zx_oq(i)
1599         zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &         zx_nq(i) = beta(i) * zx_k1(i) * (-1. * zx_dq_s_dt(i)) &
1600              &                              / zx_oq(i)              / zx_oq(i)
1601    
1602         ! H         ! H
1603         zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)         zx_oh(i) = 1. - (zx_k1(i) * petBcoef(i) * dtime)
# Line 2366  CONTAINS Line 1613  CONTAINS
1613         snow = MAX(0.0, snow)         snow = MAX(0.0, snow)
1614      end where      end where
1615    
1616      !  bil_eau_s = bil_eau_s + (precip_rain * dtime) - (evap - snow_evap) * dtime      ! bil_eau_s = bil_eau_s + (precip_rain * dtime) - (evap - snow_evap) * dtime
1617      bil_eau_s = (precip_rain * dtime) - (evap - snow_evap) * dtime      bil_eau_s = (precip_rain * dtime) - (evap - snow_evap) * dtime
1618    
1619      !  
1620      ! Y'a-t-il fonte de neige?      ! Y'a-t-il fonte de neige?
1621      !  
1622      ffonte=0.      ffonte=0.
1623      do i = 1, knon      do i = 1, knon
1624         neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &         neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
1625              & .AND. tsurf_new(i) >= RTT)              .AND. tsurf_new(i) >= RTT)
1626         if (neige_fond) then         if (neige_fond) then
1627            fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i))            fq_fonte = MIN( MAX((tsurf_new(i)-RTT )/chasno, 0.0), snow(i))
1628            ffonte(i) = fq_fonte * RLMLT/dtime            ffonte(i) = fq_fonte * RLMLT/dtime
1629            snow(i) = max(0., snow(i) - fq_fonte)            snow(i) = max(0., snow(i) - fq_fonte)
1630            bil_eau_s(i) = bil_eau_s(i) + fq_fonte            bil_eau_s(i) = bil_eau_s(i) + fq_fonte
1631            tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno              tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno
1632            !IM cf JLD OK                !IM cf JLD OK
1633            !IM cf JLD/ GKtest fonte aussi pour la glace            !IM cf JLD/ GKtest fonte aussi pour la glace
1634            IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN            IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
1635               fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,0.0)               fq_fonte = MAX((tsurf_new(i)-RTT )/chaice, 0.0)
1636               ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime               ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime
1637               bil_eau_s(i) = bil_eau_s(i) + fq_fonte               bil_eau_s(i) = bil_eau_s(i) + fq_fonte
1638               tsurf_new(i) = RTT               tsurf_new(i) = RTT
1639            ENDIF            ENDIF
1640            d_ts(i) = tsurf_new(i) - tsurf(i)            d_ts(i) = tsurf_new(i) - tsurf(i)
1641         endif         endif
1642         !  
1643         !   s'il y a une hauteur trop importante de neige, elle s'coule         ! s'il y a une hauteur trop importante de neige, elle s'coule
1644         fqcalving(i) = max(0., snow(i) - snow_max)/dtime         fqcalving(i) = max(0., snow(i) - snow_max)/dtime
1645         snow(i)=min(snow(i),snow_max)         snow(i)=min(snow(i), snow_max)
1646         !  
1647         IF (nisurf == is_ter) then         IF (nisurf == is_ter) then
1648            qsol(i) = qsol(i) + bil_eau_s(i)            qsol(i) = qsol(i) + bil_eau_s(i)
1649            run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)            run_off(i) = run_off(i) + MAX(qsol(i) - max_eau_sol, 0.0)
1650            qsol(i) = MIN(qsol(i), max_eau_sol)            qsol(i) = MIN(qsol(i), max_eau_sol)
1651         else if (nisurf == is_lic) then         else if (nisurf == is_lic) then
1652            run_off_lic(i) = (coeff_rel *  fqcalving(i)) + &            run_off_lic(i) = (coeff_rel * fqcalving(i)) + &
1653                 &                        (1. - coeff_rel) * run_off_lic_0(i)                 (1. - coeff_rel) * run_off_lic_0(i)
1654            run_off_lic_0(i) = run_off_lic(i)            run_off_lic_0(i) = run_off_lic(i)
1655            run_off_lic(i) = run_off_lic(i) + bil_eau_s(i)/dtime            run_off_lic(i) = run_off_lic(i) + bil_eau_s(i)/dtime
1656         endif         endif
Legend:
Removed from v.3  
changed lines
  Added in v.53
  ViewVC Help
Powered by ViewVC 1.1.21