/[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 14 by guez, Mon Jul 28 14:48:09 2008 UTC revision 53 by guez, Fri Oct 7 13:11:58 2011 UTC
# Line 2  MODULE interface_surf Line 2  MODULE interface_surf
2    
3    ! From phylmd/interface_surf.F90, version 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 flux de      ! en général (sols continentaux, océans, glaces) pour les flux de
# Line 62  CONTAINS Line 49  CONTAINS
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      ! input:      ! input:
57      !   klon         nombre total de points de grille      ! klon nombre total de points de grille
58      !   iim, jjm     nbres de pts de grille      ! iim, jjm nbres de pts de grille
59      !   dtime        pas de temps de la physique (en s)      ! dtime pas de temps de la physique (en s)
60      !   date0        jour initial      ! date0 jour initial
61      !   jour         jour dans l'annee en cours,      ! jour jour dans l'annee en cours,
62      !   rmu0         cosinus de l'angle solaire zenithal      ! rmu0 cosinus de l'angle solaire zenithal
63      !   nexca        pas de temps couplage      ! nexca pas de temps couplage
64      !   nisurf       index de la surface a traiter (1 = sol continental)      ! nisurf index de la surface a traiter (1 = sol continental)
65      !   knon         nombre de points de la surface a traiter      ! knon nombre de points de la surface a traiter
66      !   knindex      index des points de la surface a traiter      ! knindex index des points de la surface a traiter
67      !   pctsrf       tableau des pourcentages de surface de chaque maille      ! pctsrf tableau des pourcentages de surface de chaque maille
68      !   rlon         longitudes      ! rlon longitudes
69      !   rlat         latitudes      ! rlat latitudes
70      !   cufi, cvfi    resolution des mailles en x et y (m)      ! cufi, cvfi resolution des mailles en x et y (m)
71      !   debut        logical: 1er appel a la physique      ! debut logical: 1er appel a la physique
72      !   lafin        logical: dernier appel a la physique      ! lafin logical: dernier appel a la physique
73      !   ok_veget     logical: appel ou non au schema de surface continental      ! ok_veget logical: appel ou non au schema de surface continental
74      !                  (si false calcul simplifie des fluxs sur les continents)      ! (si false calcul simplifie des fluxs sur les continents)
75      !   zlev         hauteur de la premiere couche      ! zlev hauteur de la premiere couche
76      !   u1_lay       vitesse u 1ere couche      ! u1_lay vitesse u 1ere couche
77      !   v1_lay       vitesse v 1ere couche      ! v1_lay vitesse v 1ere couche
78      !   temp_air     temperature de l'air 1ere couche      ! temp_air temperature de l'air 1ere couche
79      !   spechum      humidite specifique 1ere couche      ! spechum humidite specifique 1ere couche
80      !   epot_air     temp potentielle de l'air      ! epot_air temp potentielle de l'air
81      !   ccanopy      concentration CO2 canopee      ! ccanopy concentration CO2 canopee
82      !   tq_cdrag     cdrag      ! tq_cdrag cdrag
83      !   petAcoef     coeff. A de la resolution de la CL pour t      ! petAcoef coeff. A de la resolution de la CL pour t
84      !   peqAcoef     coeff. A de la resolution de la CL pour q      ! peqAcoef coeff. A de la resolution de la CL pour q
85      !   petBcoef     coeff. B de la resolution de la CL pour t      ! petBcoef coeff. B de la resolution de la CL pour t
86      !   peqBcoef     coeff. B de la resolution de la CL pour q      ! peqBcoef coeff. B de la resolution de la CL pour q
87      !   precip_rain  precipitation liquide      ! precip_rain precipitation liquide
88      !   precip_snow  precipitation solide      ! precip_snow precipitation solide
89      !   sollw        flux IR net a la surface      ! sollw flux IR net a la surface
90      !   sollwdown    flux IR descendant a la surface      ! sollwdown flux IR descendant a la surface
91      !   swnet        flux solaire net      ! swnet flux solaire net
92      !   swdown       flux solaire entrant a la surface      ! swdown flux solaire entrant a la surface
93      !   albedo       albedo de la surface      ! albedo albedo de la surface
94      !   tsurf        temperature de surface      ! tsurf temperature de surface
95      !   tslab        temperature slab ocean      ! tslab temperature slab ocean
96      !   pctsrf_slab  pourcentages (0-1) des sous-surfaces dans le slab      ! pctsrf_slab pourcentages (0-1) des sous-surfaces dans le slab
97      !   tmp_pctsrf_slab = pctsrf_slab      ! tmp_pctsrf_slab = pctsrf_slab
98      !   p1lay        pression 1er niveau (milieu de couche)      ! p1lay pression 1er niveau (milieu de couche)
99      !   ps           pression au sol      ! ps pression au sol
100      !   radsol       rayonnement net aus sol (LW + SW)      ! radsol rayonnement net aus sol (LW + SW)
101      !   ocean        type d'ocean utilise ("force" ou "slab" mais pas "couple")      ! ocean type d'ocean utilise ("force" ou "slab" mais pas "couple")
102      !   fder         derivee des flux (pour le couplage)      ! fder derivee des flux (pour le couplage)
103      !   taux, tauy   tension de vents      ! taux, tauy tension de vents
104      !   windsp       module du vent a 10m      ! windsp module du vent a 10m
105      !   rugos        rugosite      ! rugos rugosite
106      !   zmasq        masque terre/ocean      ! zmasq masque terre/ocean
107      !   rugoro       rugosite orographique      ! rugoro rugosite orographique
108      !   run_off_lic_0 runoff glacier du pas de temps precedent      ! run_off_lic_0 runoff glacier du pas de temps precedent
109      integer, intent(IN) :: itime     !  numero du pas de temps      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
# Line 148  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=*), intent(IN):: 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, intent(in):: 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:      ! output:
156      !   evap         evaporation totale      ! evap evaporation totale
157      !   fluxsens     flux de chaleur sensible      ! fluxsens flux de chaleur sensible
158      !   fluxlat      flux de chaleur latente      ! fluxlat flux de chaleur latente
159      !   tsol_rad          ! tsol_rad
160      !   tsurf_new    temperature au sol      ! tsurf_new temperature au sol
161      !   alb_new      albedo      ! alb_new albedo
162      !   emis_new     emissivite      ! emis_new emissivite
163      !   z0_new       surface roughness      ! z0_new surface roughness
164      !   pctsrf_new   nouvelle repartition des surfaces      ! 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
# Line 185  CONTAINS Line 172  CONTAINS
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 203  CONTAINS Line 190  CONTAINS
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      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      :: 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
# Line 275  CONTAINS Line 262  CONTAINS
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)
# Line 285  CONTAINS Line 272  CONTAINS
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'
# Line 351  CONTAINS Line 338  CONTAINS
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    
# Line 400  CONTAINS Line 387  CONTAINS
387               CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, &               CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, &
388                    soilflux)                    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
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)
# Line 449  CONTAINS Line 436  CONTAINS
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))
# Line 463  CONTAINS Line 450  CONTAINS
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
455    
456         !IM: flux ocean-atmosphere utile pour le "slab" ocean         !IM: flux ocean-atmosphere utile pour le "slab" ocean
# Line 478  CONTAINS Line 465  CONTAINS
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)
467    
468         if (ocean == 'slab  ') then         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
# Line 494  CONTAINS Line 481  CONTAINS
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
# Line 509  CONTAINS Line 496  CONTAINS
496         ! Surface "glace de mer" appel a l'interface avec l'ocean         ! Surface "glace de mer" appel a l'interface avec l'ocean
497    
498    
499         if (ocean == 'slab  ') then         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
# Line 526  CONTAINS Line 513  CONTAINS
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 536  CONTAINS Line 523  CONTAINS
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
# Line 548  CONTAINS Line 535  CONTAINS
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
# Line 557  CONTAINS Line 544  CONTAINS
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 573  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         CALL fonte_neige( klon, knon, nisurf, dtime, &         CALL fonte_neige( klon, knon, nisurf, dtime, &
593              &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &              tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
594              &   precip_rain, precip_snow, snow, qsol,  &              precip_rain, precip_snow, snow, qsol, &
595              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
596              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              petAcoef, peqAcoef, petBcoef, peqBcoef, &
597              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &              tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
598              &   fqcalving, ffonte, run_off_lic_0)              fqcalving, ffonte, run_off_lic_0)
599    
600         !     calcul albedo         ! calcul albedo
601    
602         CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)           CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)
603         WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.         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)))         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) + &         alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
606              0.6 * (1.0-zfra(1:knon))              0.6 * (1.0-zfra(1:knon))
607    
608         fder_prev = fder             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))
# Line 626  CONTAINS Line 613  CONTAINS
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
618    
619    
# Line 654  CONTAINS Line 641  CONTAINS
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 663  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
# Line 718  CONTAINS Line 705  CONTAINS
705    
706    !************************    !************************
707    
708    SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &    SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour,  &
709         & radsol, fluxo, fluxg, pctsrf, &         radsol, fluxo, fluxg, pctsrf,  &
710         & tslab, seaice, pctsrf_slab)         tslab, seaice, pctsrf_slab)
711    
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
# Line 729  CONTAINS Line 716  CONTAINS
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 760  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 :
# Line 769  CONTAINS Line 756  CONTAINS
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
# Line 798  CONTAINS Line 785  CONTAINS
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*, '************************'
# Line 824  CONTAINS Line 811  CONTAINS
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
# Line 863  CONTAINS Line 850  CONTAINS
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    
# Line 875  CONTAINS Line 862  CONTAINS
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    
# Line 886  CONTAINS Line 873  CONTAINS
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
# Line 902  CONTAINS Line 889  CONTAINS
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      ! Cette routine sert d'interface entre le modele atmospherique et
898      ! un fichier de conditions aux limites      ! un fichier de conditions aux limites
# Line 916  CONTAINS Line 903  CONTAINS
903      use indicesol      use indicesol
904    
905      integer, intent(IN) :: itime ! numero du pas de temps courant      integer, intent(IN) :: itime ! numero du pas de temps courant
906      real   , intent(IN) :: dtime ! pas de temps de la physique (en s)      real , intent(IN) :: dtime ! pas de temps de la physique (en s)
907      integer, intent(IN) :: jour ! jour a lire dans l'annee      integer, intent(IN) :: jour ! jour a lire dans l'annee
908      integer, intent(IN) :: nisurf ! index de la surface a traiter (1 = sol continental)      integer, intent(IN) :: nisurf ! index de la surface a traiter (1 = sol continental)
909      integer, intent(IN) :: knon ! nombre de points dans le domaine a traiter      integer, intent(IN) :: knon ! nombre de points dans le domaine a traiter
# Line 926  CONTAINS Line 913  CONTAINS
913    
914      ! Parametres de sortie      ! Parametres de sortie
915      ! output:      ! output:
916      !   lmt_sst      SST lues dans le fichier de CL      ! lmt_sst SST lues dans le fichier de CL
917      !   pctsrf_new   sous-maille fractionnelle      ! 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
# Line 950  CONTAINS Line 937  CONTAINS
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      ! --------------------------------------------------      ! --------------------------------------------------
# Line 1042  CONTAINS Line 1029  CONTAINS
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
# Line 1064  CONTAINS Line 1051  CONTAINS
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)
# Line 1107  CONTAINS Line 1094  CONTAINS
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.
# Line 1121  CONTAINS Line 1108  CONTAINS
1108    
1109      ! Parametres d'entree      ! Parametres d'entree
1110      ! input:      ! input:
1111      !   itime        numero du pas de temps courant      ! itime numero du pas de temps courant
1112      !   dtime        pas de temps de la physique (en s)      ! dtime pas de temps de la physique (en s)
1113      !   jour         jour a lire dans l'annee      ! jour jour a lire dans l'annee
1114      !   nisurf       index de la surface a traiter (1 = sol continental)      ! nisurf index de la surface a traiter (1 = sol continental)
1115      !   knon         nombre de points dans le domaine a traiter      ! knon nombre de points dans le domaine a traiter
1116      !   knindex      index des points de la surface a traiter      ! knindex index des points de la surface a traiter
1117      !   klon         taille de la grille      ! klon taille de la grille
1118      !   debut        logical: 1er appel a la physique (initialisation)      ! 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 1140  CONTAINS Line 1127  CONTAINS
1127    
1128      ! Parametres de sortie      ! Parametres de sortie
1129      ! output:      ! output:
1130      !   lmt_sst      SST lues dans le fichier de CL      ! lmt_sst SST lues dans le fichier de CL
1131      !   lmt_alb      Albedo lu      ! lmt_alb Albedo lu
1132      !   lmt_rug      longueur de rugosité lue      ! lmt_rug longueur de rugosité lue
1133      !   pctsrf_new   sous-maille fractionnelle      ! 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      !------------------------------------------------------------      !------------------------------------------------------------
# Line 1185  CONTAINS Line 1172  CONTAINS
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 &
# Line 1201  CONTAINS Line 1188  CONTAINS
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'
# Line 1214  CONTAINS Line 1201  CONTAINS
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'
# Line 1228  CONTAINS Line 1215  CONTAINS
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
# Line 1238  CONTAINS Line 1225  CONTAINS
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
# Line 1251  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)
# Line 1264  CONTAINS Line 1251  CONTAINS
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 1321  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    
# Line 1353  CONTAINS Line 1340  CONTAINS
1340    
1341      ! Traitement neige et humidite du sol      ! Traitement neige et humidite du sol
1342    
 !!$  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  
   
   
1343      ! Initialisation      ! Initialisation
1344    
1345      evap = 0.      evap = 0.
# Line 1388  CONTAINS Line 1361  CONTAINS
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    
# Line 1422  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 1435  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 1469  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 1542  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 1551  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 1563  CONTAINS Line 1536  CONTAINS
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    
# Line 1587  CONTAINS Line 1560  CONTAINS
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 ===
# Line 1620  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 1640  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    
# Line 1649  CONTAINS Line 1622  CONTAINS
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)
# Line 1667  CONTAINS Line 1640  CONTAINS
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    
# Line 1676  CONTAINS Line 1649  CONTAINS
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.14  
changed lines
  Added in v.53
  ViewVC Help
Powered by ViewVC 1.1.21