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

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

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

revision 12 by guez, Mon Jul 21 16:05:07 2008 UTC revision 14 by guez, Mon Jul 28 14:48:09 2008 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 modèle
6    ! atmosphérique et les modèles de surface (sols continentaux,    ! atmosphérique et les modèles de surface (sols continentaux,
# Line 18  MODULE interface_surf Line 18  MODULE interface_surf
18    PUBLIC :: interfsurf_hq    PUBLIC :: interfsurf_hq
19    
20    ! run_off      ruissellement total    ! run_off      ruissellement total
21    REAL, ALLOCATABLE, DIMENSION(:),SAVE    :: run_off, run_off_lic    REAL, ALLOCATABLE, DIMENSION(:), SAVE    :: run_off, run_off_lic
22    real, allocatable, dimension(:),save    :: coastalflow, riverflow    real, allocatable, dimension(:), save    :: coastalflow, riverflow
23    
24    REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: tmp_rriv, tmp_rcoa,tmp_rlic    REAL, ALLOCATABLE, DIMENSION(:, :), SAVE :: tmp_rriv, tmp_rcoa, tmp_rlic
25    !! pour simuler la fonte des glaciers antarctiques    !! pour simuler la fonte des glaciers antarctiques
26    REAL, ALLOCATABLE, DIMENSION(:,:), SAVE :: coeff_iceberg    REAL, ALLOCATABLE, DIMENSION(:, :), SAVE :: coeff_iceberg
27    real, save                              :: surf_maille    real, save                              :: surf_maille
28    real, save                              :: cte_flux_iceberg = 6.3e7    real, save                              :: cte_flux_iceberg = 6.3e7
29    integer, save                           :: num_antarctic = 1    integer, save                           :: num_antarctic = 1
# Line 33  CONTAINS Line 33  CONTAINS
33    
34    SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, &    SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, &
35         & klon, iim, jjm, nisurf, knon, knindex, pctsrf, &         & klon, iim, jjm, nisurf, knon, knindex, pctsrf, &
36         & rlon, rlat, cufi, cvfi,&         & rlon, rlat, cufi, cvfi, &
37         & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, qsol,&         & debut, lafin, ok_veget, soil_model, nsoilmx, tsoil, qsol, &
38         & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &         & zlev,  u1_lay, v1_lay, temp_air, spechum, epot_air, ccanopy, &
39         & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &         & tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
40         & precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &         & precip_rain, precip_snow, sollw, sollwdown, swnet, swdown, &
# Line 46  CONTAINS Line 46  CONTAINS
46         & ocean, npas, nexca, zmasq, &         & ocean, npas, nexca, zmasq, &
47         & evap, fluxsens, fluxlat, dflux_l, dflux_s, &                       & evap, fluxsens, fluxlat, dflux_l, dflux_s, &              
48         & tsol_rad, tsurf_new, alb_new, alblw, emis_new, &         & tsol_rad, tsurf_new, alb_new, alblw, emis_new, &
49         & z0_new, pctsrf_new, agesno,fqcalving,ffonte, run_off_lic_0,&         & z0_new, pctsrf_new, agesno, fqcalving, ffonte, run_off_lic_0, &
50         !IM "slab" ocean         !IM "slab" ocean
51         & flux_o, flux_g, tslab, seaice)         & flux_o, flux_g, tslab, seaice)
52    
53      ! Cette routine sert d'aiguillage entre l'atmosphère et la surface      ! Cette routine sert d'aiguillage entre l'atmosphère et la surface
54      ! 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
55      ! chaleur et d'humidité.      ! chaleur et d'humidité.
56      ! 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
57      ! atmosphérique ("clmain.F") et les routines de surface      ! atmosphérique ("clmain.F") et les routines de surface
# Line 59  CONTAINS Line 59  CONTAINS
59    
60      ! L.Fairhead 02/2000      ! L.Fairhead 02/2000
61    
62        use abort_gcm_m, only: abort_gcm
63        use gath_cpl, only: gath2cpl
64        use indicesol
65        use YOMCST
66        use albsno_m, only: albsno
67    
68        ! Parametres d'entree
69      ! input:      ! input:
70      !   klon         nombre total de points de grille      !   klon         nombre total de points de grille
71      !   iim, jjm     nbres de pts de grille      !   iim, jjm     nbres de pts de grille
# Line 73  CONTAINS Line 80  CONTAINS
80      !   pctsrf       tableau des pourcentages de surface de chaque maille      !   pctsrf       tableau des pourcentages de surface de chaque maille
81      !   rlon         longitudes      !   rlon         longitudes
82      !   rlat         latitudes      !   rlat         latitudes
83      !   cufi,cvfi    resolution des mailles en x et y (m)      !   cufi, cvfi    resolution des mailles en x et y (m)
84      !   debut        logical: 1er appel a la physique      !   debut        logical: 1er appel a la physique
85      !   lafin        logical: dernier appel a la physique      !   lafin        logical: dernier appel a la physique
86      !   ok_veget     logical: appel ou non au schema de surface continental      !   ok_veget     logical: appel ou non au schema de surface continental
# Line 112  CONTAINS Line 119  CONTAINS
119      !   zmasq        masque terre/ocean      !   zmasq        masque terre/ocean
120      !   rugoro       rugosite orographique      !   rugoro       rugosite orographique
121      !   run_off_lic_0 runoff glacier du pas de temps precedent      !   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  
   
     use abort_gcm_m, only: abort_gcm  
     use gath_cpl, only: gath2cpl  
     use indicesol  
     use YOMCST  
     use albsno_m, only: albsno  
   
     ! Parametres d'entree  
122      integer, intent(IN) :: itime     !  numero du pas de temps      integer, intent(IN) :: itime     !  numero du pas de temps
123      integer, intent(IN) :: iim, jjm      integer, intent(IN) :: iim, jjm
124      integer, intent(IN) :: klon      integer, intent(IN) :: klon
# Line 141  CONTAINS Line 129  CONTAINS
129      integer, intent(IN) :: nisurf      integer, intent(IN) :: nisurf
130      integer, intent(IN) :: knon      integer, intent(IN) :: knon
131      integer, dimension(klon), intent(in) :: knindex      integer, dimension(klon), intent(in) :: knindex
132      real, dimension(klon,nbsrf), intent(IN) :: pctsrf      real, dimension(klon, nbsrf), intent(IN) :: pctsrf
133      logical, intent(IN) :: debut, lafin, ok_veget      logical, intent(IN) :: debut, lafin, ok_veget
134      real, dimension(klon), intent(IN) :: rlon, rlat      real, dimension(klon), intent(IN) :: rlon, rlat
135      real, dimension(klon), intent(IN) :: cufi, cvfi      real, dimension(klon), intent(IN) :: cufi, cvfi
# Line 161  CONTAINS Line 149  CONTAINS
149      real, allocatable, dimension(:), save :: tmp_tslab      real, allocatable, dimension(:), save :: tmp_tslab
150      real, dimension(klon), intent(OUT) :: flux_o, flux_g      real, dimension(klon), intent(OUT) :: flux_o, flux_g
151      real, dimension(klon), intent(INOUT)      :: seaice ! glace de mer (kg/m2)      real, dimension(klon), intent(INOUT)      :: seaice ! glace de mer (kg/m2)
152      REAL, DIMENSION(klon), INTENT(INOUT) :: radsol,fder      REAL, DIMENSION(klon), INTENT(INOUT) :: radsol, fder
153      real, dimension(klon), intent(IN) :: zmasq      real, dimension(klon), intent(IN) :: zmasq
154      real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro      real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro
155      real, dimension(klon), intent(IN) :: windsp      real, dimension(klon), intent(IN) :: windsp
# Line 175  CONTAINS Line 163  CONTAINS
163      REAL, dimension(klon), intent(INOUT) :: qsol      REAL, dimension(klon), intent(INOUT) :: qsol
164      REAL, dimension(klon)          :: soilcap      REAL, dimension(klon)          :: soilcap
165      REAL, dimension(klon)          :: soilflux      REAL, dimension(klon)          :: soilflux
166    
167      ! Parametres de sortie      ! Parametres de sortie
168        ! output:
169        !   evap         evaporation totale
170        !   fluxsens     flux de chaleur sensible
171        !   fluxlat      flux de chaleur latente
172        !   tsol_rad    
173        !   tsurf_new    temperature au sol
174        !   alb_new      albedo
175        !   emis_new     emissivite
176        !   z0_new       surface roughness
177        !   pctsrf_new   nouvelle repartition des surfaces
178      real, dimension(klon), intent(OUT):: fluxsens, fluxlat      real, dimension(klon), intent(OUT):: fluxsens, fluxlat
179      real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new      real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new
180      real, dimension(klon), intent(OUT):: alblw      real, dimension(klon), intent(OUT):: alblw
181      real, dimension(klon), intent(OUT):: emis_new, z0_new      real, dimension(klon), intent(OUT):: emis_new, z0_new
182      real, dimension(klon), intent(OUT):: dflux_l, dflux_s      real, dimension(klon), intent(OUT):: dflux_l, dflux_s
183      real, dimension(klon,nbsrf), intent(OUT) :: pctsrf_new      real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_new
184      real, dimension(klon), intent(INOUT):: agesno      real, dimension(klon), intent(INOUT):: agesno
185      real, dimension(klon), intent(INOUT):: run_off_lic_0      real, dimension(klon), intent(INOUT):: run_off_lic_0
186    
# Line 198  CONTAINS Line 197  CONTAINS
197      integer i      integer i
198      real, allocatable, dimension(:), save :: tmp_flux_o, tmp_flux_g      real, allocatable, dimension(:), save :: tmp_flux_o, tmp_flux_g
199      real, allocatable, dimension(:), save :: tmp_radsol      real, allocatable, dimension(:), save :: tmp_radsol
200      real, allocatable, dimension(:,:), save :: tmp_pctsrf_slab      real, allocatable, dimension(:, :), save :: tmp_pctsrf_slab
201      real, allocatable, dimension(:), save :: tmp_seaice      real, allocatable, dimension(:), save :: tmp_seaice
202    
203      ! Local      ! Local
204      character (len = 20),save :: modname = 'interfsurf_hq'      character (len = 20), save :: modname = 'interfsurf_hq'
205      character (len = 80) :: abort_message      character (len = 80) :: abort_message
206      logical, save        :: first_call = .true.      logical, save        :: first_call = .true.
207      integer, save        :: error      integer, save        :: error
208      integer              :: ii      integer              :: ii
209      logical,save              :: check = .false.      logical, save              :: check = .false.
210      real, dimension(klon):: cal, beta, dif_grnd, capsol      real, dimension(klon):: cal, beta, dif_grnd, capsol
 !!$PB  real, parameter      :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.  
211      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.
212      real, parameter      :: calsno=1./(2.3867e+06*.15)      real, parameter      :: calsno=1./(2.3867e+06*.15)
213      real, dimension(klon):: tsurf_temp      real, dimension(klon):: tsurf_temp
214      real, dimension(klon):: alb_neig, alb_eau      real, dimension(klon):: alb_neig, alb_eau
215      real, DIMENSION(klon):: zfra      real, DIMENSION(klon):: zfra
216      logical              :: cumul = .false.      logical              :: cumul = .false.
217      INTEGER,dimension(1) :: iloc      INTEGER, dimension(1) :: iloc
218      real, dimension(klon):: fder_prev      real, dimension(klon):: fder_prev
219      REAL, dimension(klon) :: bidule      REAL, dimension(klon) :: bidule
220    
221      !-------------------------------------------------------------      !-------------------------------------------------------------
222    
223      if (check) write(*,*) 'Entree ', modname      if (check) write(*, *) 'Entree ', modname
224    
225      ! On doit commencer par appeler les schemas de surfaces continentales      ! On doit commencer par appeler les schemas de surfaces continentales
226      ! 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 228  CONTAINS
228      if (first_call) then      if (first_call) then
229         call conf_interface(tau_calv)         call conf_interface(tau_calv)
230         if (nisurf /= is_ter .and. klon > 1) then         if (nisurf /= is_ter .and. klon > 1) then
231            write(*,*)' *** Warning ***'            write(*, *)' *** Warning ***'
232            write(*,*)' nisurf = ',nisurf,' /= is_ter = ',is_ter            write(*, *)' nisurf = ', nisurf, ' /= is_ter = ', is_ter
233            write(*,*)'or on doit commencer par les surfaces continentales'            write(*, *)'or on doit commencer par les surfaces continentales'
234            abort_message='voir ci-dessus'            abort_message='voir ci-dessus'
235            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
236         endif         endif
237         if (ocean /= 'slab' .and. ocean /= 'force') then         if (ocean /= 'slab' .and. ocean /= 'force') then
238            write(*,*)' *** Warning ***'            write(*, *)' *** Warning ***'
239            write(*,*)'Option couplage pour l''ocean = ', ocean            write(*, *)'Option couplage pour l''ocean = ', ocean
240            abort_message='option pour l''ocean non valable'            abort_message='option pour l''ocean non valable'
241            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
242         endif         endif
243         if ( is_oce > is_sic ) then         if ( is_oce > is_sic ) then
244            write(*,*)' *** Warning ***'            write(*, *)' *** Warning ***'
245            write(*,*)' Pour des raisons de sequencement dans le code'            write(*, *)' Pour des raisons de sequencement dans le code'
246            write(*,*)' l''ocean doit etre traite avant la banquise'            write(*, *)' l''ocean doit etre traite avant la banquise'
247            write(*,*)' or is_oce = ',is_oce, '> is_sic = ',is_sic            write(*, *)' or is_oce = ', is_oce, '> is_sic = ', is_sic
248            abort_message='voir ci-dessus'            abort_message='voir ci-dessus'
249            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
250         endif         endif
251      endif      endif
252      first_call = .false.      first_call = .false.
253    
254      ! Initialisations diverses      ! Initialisations diverses
255      !  
256      ffonte(1:knon)=0.      ffonte(1:knon)=0.
257      fqcalving(1:knon)=0.      fqcalving(1:knon)=0.
258    
# Line 266  CONTAINS Line 264  CONTAINS
264      !IM: "slab" ocean; initialisations      !IM: "slab" ocean; initialisations
265      flux_o = 0.      flux_o = 0.
266      flux_g = 0.      flux_g = 0.
267      !  
268      if (.not. allocated(tmp_flux_o)) then      if (.not. allocated(tmp_flux_o)) then
269         allocate(tmp_flux_o(klon), stat = error)         allocate(tmp_flux_o(klon), stat = error)
270         DO i=1, knon         DO i=1, knon
# Line 274  CONTAINS Line 272  CONTAINS
272         ENDDO         ENDDO
273         if (error /= 0) then         if (error /= 0) then
274            abort_message='Pb allocation tmp_flux_o'            abort_message='Pb allocation tmp_flux_o'
275            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
276         endif         endif
277      endif      endif
278      if (.not. allocated(tmp_flux_g)) then        if (.not. allocated(tmp_flux_g)) then  
# Line 284  CONTAINS Line 282  CONTAINS
282         ENDDO         ENDDO
283         if (error /= 0) then         if (error /= 0) then
284            abort_message='Pb allocation tmp_flux_g'            abort_message='Pb allocation tmp_flux_g'
285            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
286         endif         endif
287      endif      endif
288      if (.not. allocated(tmp_radsol)) then        if (.not. allocated(tmp_radsol)) then  
289         allocate(tmp_radsol(klon), stat = error)         allocate(tmp_radsol(klon), stat = error)
290         if (error /= 0) then         if (error /= 0) then
291            abort_message='Pb allocation tmp_radsol'            abort_message='Pb allocation tmp_radsol'
292            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
293         endif         endif
294      endif      endif
295      DO i=1, knon      DO i=1, knon
296         tmp_radsol(knindex(i))=radsol(i)         tmp_radsol(knindex(i))=radsol(i)
297      ENDDO      ENDDO
298      if (.not. allocated(tmp_pctsrf_slab)) then      if (.not. allocated(tmp_pctsrf_slab)) then
299         allocate(tmp_pctsrf_slab(klon,nbsrf), stat = error)         allocate(tmp_pctsrf_slab(klon, nbsrf), stat = error)
300         if (error /= 0) then         if (error /= 0) then
301            abort_message='Pb allocation tmp_pctsrf_slab'            abort_message='Pb allocation tmp_pctsrf_slab'
302            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
303         endif         endif
304         DO i=1, klon         DO i=1, klon
305            tmp_pctsrf_slab(i,1:nbsrf)=pctsrf(i,1:nbsrf)            tmp_pctsrf_slab(i, 1:nbsrf)=pctsrf(i, 1:nbsrf)
306         ENDDO         ENDDO
307      endif      endif
308      !  
309      if (.not. allocated(tmp_seaice)) then      if (.not. allocated(tmp_seaice)) then
310         allocate(tmp_seaice(klon), stat = error)         allocate(tmp_seaice(klon), stat = error)
311         if (error /= 0) then         if (error /= 0) then
312            abort_message='Pb allocation tmp_seaice'            abort_message='Pb allocation tmp_seaice'
313            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
314         endif         endif
315         DO i=1, klon         DO i=1, klon
316            tmp_seaice(i)=seaice(i)            tmp_seaice(i)=seaice(i)
317         ENDDO         ENDDO
318      endif      endif
319      !  
320      if (.not. allocated(tmp_tslab)) then      if (.not. allocated(tmp_tslab)) then
321         allocate(tmp_tslab(klon), stat = error)         allocate(tmp_tslab(klon), stat = error)
322         if (error /= 0) then         if (error /= 0) then
323            abort_message='Pb allocation tmp_tslab'            abort_message='Pb allocation tmp_tslab'
324            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
325         endif         endif
326      endif      endif
327      DO i=1, klon      DO i=1, klon
328         tmp_tslab(i)=tslab(i)         tmp_tslab(i)=tslab(i)
329      ENDDO      ENDDO
330      !  
331      ! Aiguillage vers les differents schemas de surface      ! Aiguillage vers les differents schemas de surface
332    
333      if (nisurf == is_ter) then      if (nisurf == is_ter) then
334         !  
335         ! Surface "terre" appel a l'interface avec les sols continentaux         ! Surface "terre" appel a l'interface avec les sols continentaux
336         !  
337         ! allocation du run-off         ! allocation du run-off
338         if (.not. allocated(coastalflow)) then         if (.not. allocated(coastalflow)) then
339            allocate(coastalflow(knon), stat = error)            allocate(coastalflow(knon), stat = error)
340            if (error /= 0) then            if (error /= 0) then
341               abort_message='Pb allocation coastalflow'               abort_message='Pb allocation coastalflow'
342               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
343            endif            endif
344            allocate(riverflow(knon), stat = error)            allocate(riverflow(knon), stat = error)
345            if (error /= 0) then            if (error /= 0) then
346               abort_message='Pb allocation riverflow'               abort_message='Pb allocation riverflow'
347               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
348            endif            endif
349            allocate(run_off(knon), stat = error)            allocate(run_off(knon), stat = error)
350            if (error /= 0) then            if (error /= 0) then
351               abort_message='Pb allocation run_off'               abort_message='Pb allocation run_off'
352               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
353            endif            endif
354            !cym                  !cym      
355            run_off=0.0            run_off=0.0
356            !cym            !cym
357    
358  !!$PB  !!$PB
359            ALLOCATE (tmp_rriv(iim,jjm+1), stat=error)            ALLOCATE (tmp_rriv(iim, jjm+1), stat=error)
360            if (error /= 0) then            if (error /= 0) then
361               abort_message='Pb allocation tmp_rriv'               abort_message='Pb allocation tmp_rriv'
362               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
363            endif            endif
364            ALLOCATE (tmp_rcoa(iim,jjm+1), stat=error)            ALLOCATE (tmp_rcoa(iim, jjm+1), stat=error)
365            if (error /= 0) then            if (error /= 0) then
366               abort_message='Pb allocation tmp_rcoa'               abort_message='Pb allocation tmp_rcoa'
367               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
368            endif            endif
369            ALLOCATE (tmp_rlic(iim,jjm+1), stat=error)            ALLOCATE (tmp_rlic(iim, jjm+1), stat=error)
370            if (error /= 0) then            if (error /= 0) then
371               abort_message='Pb allocation tmp_rlic'               abort_message='Pb allocation tmp_rlic'
372               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
373            endif            endif
374            tmp_rriv = 0.0            tmp_rriv = 0.0
375            tmp_rcoa = 0.0            tmp_rcoa = 0.0
# Line 379  CONTAINS Line 377  CONTAINS
377    
378  !!$  !!$
379         else if (size(coastalflow) /= knon) then         else if (size(coastalflow) /= knon) then
380            write(*,*)'Bizarre, le nombre de points continentaux'            write(*, *)'Bizarre, le nombre de points continentaux'
381            write(*,*)'a change entre deux appels. J''arrete ...'            write(*, *)'a change entre deux appels. J''arrete ...'
382            abort_message='voir ci-dessus'            abort_message='voir ci-dessus'
383            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
384         endif         endif
385         coastalflow = 0.         coastalflow = 0.
386         riverflow = 0.         riverflow = 0.
387         !  
388         ! Calcul age de la neige         ! Calcul age de la neige
389    
390         if (.not. ok_veget) then         if (.not. ok_veget) then
391            !            ! calcul albedo: lecture albedo fichier boundary conditions
392            ! calcul albedo: lecture albedo fichier CL puis ajout albedo neige            ! puis ajout albedo neige
393            !            call interfsur_lim(itime, dtime, jour, klon, nisurf, knon, knindex, &
394            call interfsur_lim(itime, dtime, jour, &                 debut, alb_new, z0_new)
395                 & klon, nisurf, knon, knindex, debut,  &  
                & alb_new, z0_new)  
           !    
396            ! calcul snow et qsurf, hydrol adapté            ! calcul snow et qsurf, hydrol adapté
           !  
397            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
398    
399            IF (soil_model) THEN            IF (soil_model) THEN
400               CALL soil(dtime, nisurf, knon,snow, tsurf, tsoil,soilcap, soilflux)               CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, &
401                      soilflux)
402               cal(1:knon) = RCPD / soilcap(1:knon)               cal(1:knon) = RCPD / soilcap(1:knon)
403               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
404            ELSE            ELSE
405               cal = RCPD * capsol               cal = RCPD * capsol
 !!$      cal = capsol  
406            ENDIF            ENDIF
407            CALL calcul_fluxs( klon, knon, nisurf, dtime, &            CALL calcul_fluxs( klon, knon, nisurf, dtime, &
408                 &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &                 &   tsurf, p1lay, cal, beta, tq_cdrag, ps, &
# Line 422  CONTAINS Line 417  CONTAINS
417                 &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &                 &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
418                 &   petAcoef, peqAcoef, petBcoef, peqBcoef, &                 &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
419                 &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &                 &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
420                 &   fqcalving,ffonte, run_off_lic_0)                 &   fqcalving, ffonte, run_off_lic_0)
421    
422            call albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))              call albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)
423            where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.            where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
424            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)))
425            alb_new(1 : knon)  = alb_neig(1 : knon) *zfra(1:knon) + &            alb_new(1 : knon)  = alb_neig(1 : knon) *zfra(1:knon) + &
426                 &                     alb_new(1 : knon)*(1.0-zfra(1:knon))                 alb_new(1 : knon)*(1.0-zfra(1:knon))
427            z0_new = sqrt(z0_new**2+rugoro**2)            z0_new = sqrt(z0_new**2+rugoro**2)
428            alblw(1 : knon) = alb_new(1 : knon)            alblw(1 : knon) = alb_new(1 : knon)
   
        else  
429         endif         endif
        !  
        ! Remplissage des pourcentages de surface  
        !  
        pctsrf_new(:,nisurf) = pctsrf(:,nisurf)  
430    
431           ! Remplissage des pourcentages de surface
432           pctsrf_new(:, nisurf) = pctsrf(:, nisurf)
433      else if (nisurf == is_oce) then      else if (nisurf == is_oce) then
434         ! Surface "ocean" appel a l'interface avec l'ocean         ! Surface "ocean" appel a l'interface avec l'ocean
435         !         if (ocean == 'slab') then
        if (ocean == 'slab  ') then  
436            tsurf_new = tsurf            tsurf_new = tsurf
437            pctsrf_new = tmp_pctsrf_slab            pctsrf_new = tmp_pctsrf_slab
438            !         else
439         else                              ! lecture conditions limites            ! lecture conditions limites
440            call interfoce_lim(itime, dtime, jour, &            call interfoce_lim(itime, dtime, jour, klon, nisurf, knon, knindex, &
441                 &  klon, nisurf, knon, knindex, &                 debut, tsurf_new, pctsrf_new)
                &  debut, &  
                &  tsurf_new, pctsrf_new)  
   
442         endif         endif
443    
444         tsurf_temp = tsurf_new         tsurf_temp = tsurf_new
445         cal = 0.         cal = 0.
446         beta = 1.         beta = 1.
447         dif_grnd = 0.         dif_grnd = 0.
448         alb_neig(:) = 0.         alb_neig = 0.
449         agesno(:) = 0.         agesno = 0.
450    
451         call calcul_fluxs( klon, knon, nisurf, dtime, &         call calcul_fluxs( klon, knon, nisurf, dtime, &
452              &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &              &   tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
# Line 473  CONTAINS Line 460  CONTAINS
460    
461         iloc = maxloc(fder(1:klon))         iloc = maxloc(fder(1:klon))
462         if (check.and.fder(iloc(1))> 0.) then         if (check.and.fder(iloc(1))> 0.) then
463            WRITE(*,*)'**** Debug fder****'            WRITE(*, *)'**** Debug fder****'
464            WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))            WRITE(*, *)'max fder(', iloc(1), ') = ', fder(iloc(1))
465            WRITE(*,*)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &            WRITE(*, *)'fder_prev, dflux_s, dflux_l', fder_prev(iloc(1)), &
466                 &                        dflux_s(iloc(1)), dflux_l(iloc(1))                 &                        dflux_s(iloc(1)), dflux_l(iloc(1))
467         endif         endif
 !!$  
 !!$      where(fder.gt.0.)  
 !!$        fder = 0.  
 !!$      endwhere  
468    
469         !IM: flux ocean-atmosphere utile pour le "slab" ocean         !IM: flux ocean-atmosphere utile pour le "slab" ocean
470         DO i=1, knon         DO i=1, knon
# Line 491  CONTAINS Line 474  CONTAINS
474            tmp_flux_o(knindex(i)) = flux_o(i)            tmp_flux_o(knindex(i)) = flux_o(i)
475            tmp_radsol(knindex(i))=radsol(i)            tmp_radsol(knindex(i))=radsol(i)
476         ENDDO         ENDDO
477         !  
478         ! 2eme appel a interfoce pour le cumul des champs (en particulier         ! 2eme appel a interfoce pour le cumul des champs (en particulier
479         ! fluxsens et fluxlat calcules dans calcul_fluxs)         ! fluxsens et fluxlat calcules dans calcul_fluxs)
480         !  
481         if (ocean == 'slab  ') then         if (ocean == 'slab  ') then
           !  
482            seaice=tmp_seaice            seaice=tmp_seaice
483            cumul = .true.            cumul = .true.
484            call interfoce_slab(klon, debut, itime, dtime, jour, &            call interfoce_slab(klon, debut, itime, dtime, jour, &
485                 & tmp_radsol, tmp_flux_o, tmp_flux_g, pctsrf, &                 & tmp_radsol, tmp_flux_o, tmp_flux_g, pctsrf, &
486                 & tslab, seaice, pctsrf_new)                 & tslab, seaice, pctsrf_new)
487            !  
488            tmp_pctsrf_slab=pctsrf_new            tmp_pctsrf_slab=pctsrf_new
489            DO i=1, knon            DO i=1, knon
490               tsurf_new(i)=tslab(knindex(i))               tsurf_new(i)=tslab(knindex(i))
491            ENDDO !i            ENDDO
           !  
492         endif         endif
493    
        !  
494         ! calcul albedo         ! calcul albedo
        !  
   
495         if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then         if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then
496            CALL alboc(FLOAT(jour),rlat,alb_eau)            CALL alboc(FLOAT(jour), rlat, alb_eau)
497         else  ! cycle diurne         else  ! cycle diurne
498            CALL alboc_cd(rmu0,alb_eau)            CALL alboc_cd(rmu0, alb_eau)
499         endif         endif
500         DO ii =1, knon         DO ii =1, knon
501            alb_new(ii) = alb_eau(knindex(ii))            alb_new(ii) = alb_eau(knindex(ii))
# Line 525  CONTAINS Line 503  CONTAINS
503    
504         z0_new = sqrt(rugos**2 + rugoro**2)         z0_new = sqrt(rugos**2 + rugoro**2)
505         alblw(1:knon) = alb_new(1:knon)         alblw(1:knon) = alb_new(1:knon)
   
        !  
506      else if (nisurf == is_sic) then      else if (nisurf == is_sic) then
507           if (check) write(*, *)'sea ice, nisurf = ', nisurf
508    
        if (check) write(*,*)'sea ice, nisurf = ',nisurf  
   
        !  
509         ! Surface "glace de mer" appel a l'interface avec l'ocean         ! Surface "glace de mer" appel a l'interface avec l'ocean
510         !  
511         !  
512         if (ocean == 'slab  ') then         if (ocean == 'slab  ') then
513            pctsrf_new=tmp_pctsrf_slab            pctsrf_new=tmp_pctsrf_slab
514            !  
515            DO ii = 1, knon            DO ii = 1, knon
516               tsurf_new(ii) = tsurf(ii)               tsurf_new(ii) = tsurf(ii)
517               IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then               IF (pctsrf_new(knindex(ii), nisurf) < EPSFRA) then
518                  snow(ii) = 0.0                  snow(ii) = 0.0
519                  tsurf_new(ii) = RTT - 1.8                  tsurf_new(ii) = RTT - 1.8
520                  IF (soil_model) tsoil(ii,:) = RTT -1.8                  IF (soil_model) tsoil(ii, :) = RTT -1.8
521               ENDIF               ENDIF
522            ENDDO            ENDDO
523    
524            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
525    
526            IF (soil_model) THEN            IF (soil_model) THEN
527               CALL soil(dtime, nisurf, knon,snow, tsurf_new, tsoil,soilcap, soilflux)               CALL soil(dtime, nisurf, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
528               cal(1:knon) = RCPD / soilcap(1:knon)               cal(1:knon) = RCPD / soilcap(1:knon)
529               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
530            ELSE            ELSE
# Line 560  CONTAINS Line 534  CONTAINS
534            ENDIF            ENDIF
535            tsurf_temp = tsurf_new            tsurf_temp = tsurf_new
536            beta = 1.0            beta = 1.0
537            !  
538         ELSE         ELSE
539            !                              ! lecture conditions limites            !                              ! lecture conditions limites
540            CALL interfoce_lim(itime, dtime, jour, &            CALL interfoce_lim(itime, dtime, jour, &
# Line 571  CONTAINS Line 545  CONTAINS
545            !IM cf LF            !IM cf LF
546            DO ii = 1, knon            DO ii = 1, knon
547               tsurf_new(ii) = tsurf(ii)               tsurf_new(ii) = tsurf(ii)
548               !IMbad IF (pctsrf_new(ii,nisurf) < EPSFRA) then               !IMbad IF (pctsrf_new(ii, nisurf) < EPSFRA) then
549               IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then               IF (pctsrf_new(knindex(ii), nisurf) < EPSFRA) then
550                  snow(ii) = 0.0                  snow(ii) = 0.0
551                  !IM cf LF/JLD         tsurf(ii) = RTT - 1.8                  !IM cf LF/JLD         tsurf(ii) = RTT - 1.8
552                  tsurf_new(ii) = RTT - 1.8                  tsurf_new(ii) = RTT - 1.8
553                  IF (soil_model) tsoil(ii,:) = RTT -1.8                  IF (soil_model) tsoil(ii, :) = RTT -1.8
554               endif               endif
555            enddo            enddo
556    
557            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)            CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
558    
559            IF (soil_model) THEN            IF (soil_model) THEN
560               !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)
561               CALL soil(dtime, nisurf, knon,snow, tsurf_new, tsoil,soilcap, soilflux)               CALL soil(dtime, nisurf, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
562               cal(1:knon) = RCPD / soilcap(1:knon)               cal(1:knon) = RCPD / soilcap(1:knon)
563               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)               radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
564               dif_grnd = 0.               dif_grnd = 0.
# Line 604  CONTAINS Line 578  CONTAINS
578              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
579              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
580              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
581         !  
582         !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
583         DO i = 1, knon         DO i = 1, knon
584            flux_g(i) = 0.0            flux_g(i) = 0.0
585            !  
586            !IM: faire dependre le coefficient de conduction de la glace de mer            !IM: faire dependre le coefficient de conduction de la glace de mer
587            !    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.
588            !    actuel correspond a 3m de glace de mer, cf. L.Li            !    actuel correspond a 3m de glace de mer, cf. L.Li
589            !  
590            !      IF(1.EQ.0) THEN            !      IF(1.EQ.0) THEN
591            !       IF(siceh(i).GT.0.) THEN            !       IF(siceh(i).GT.0.) THEN
592            !        new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)            !        new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)
# Line 620  CONTAINS Line 594  CONTAINS
594            !        new_dif_grnd(i) = 0.            !        new_dif_grnd(i) = 0.
595            !       ENDIF            !       ENDIF
596            !      ENDIF !(1.EQ.0) THEN            !      ENDIF !(1.EQ.0) THEN
597            !  
598            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) &
599                 &                          * dif_grnd(i) *RCPD/cal(i)                 &                          * dif_grnd(i) *RCPD/cal(i)
600            !   &                          * new_dif_grnd(i) *RCPD/cal(i)            !   &                          * new_dif_grnd(i) *RCPD/cal(i)
# Line 634  CONTAINS Line 608  CONTAINS
608              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
609              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
610              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
611              &   fqcalving,ffonte, run_off_lic_0)              &   fqcalving, ffonte, run_off_lic_0)
612    
613         !     calcul albedo         !     calcul albedo
614    
615         CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))           CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)  
616         WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.         WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
617         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)))
618         alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &         alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
619              0.6 * (1.0-zfra(1:knon))              0.6 * (1.0-zfra(1:knon))
620    
# Line 649  CONTAINS Line 623  CONTAINS
623    
624         iloc = maxloc(fder(1:klon))         iloc = maxloc(fder(1:klon))
625         if (check.and.fder(iloc(1))> 0.) then         if (check.and.fder(iloc(1))> 0.) then
626            WRITE(*,*)'**** Debug fder ****'            WRITE(*, *)'**** Debug fder ****'
627            WRITE(*,*)'max fder(',iloc(1),') = ',fder(iloc(1))            WRITE(*, *)'max fder(', iloc(1), ') = ', fder(iloc(1))
628            WRITE(*,*)'fder_prev, dflux_s, dflux_l',fder_prev(iloc(1)), &            WRITE(*, *)'fder_prev, dflux_s, dflux_l', fder_prev(iloc(1)), &
629                 &                        dflux_s(iloc(1)), dflux_l(iloc(1))                 &                        dflux_s(iloc(1)), dflux_l(iloc(1))
630         endif         endif
631    
632         !  
633         ! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean         ! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
634         !  
635         z0_new = 0.002         z0_new = 0.002
636         z0_new = SQRT(z0_new**2+rugoro**2)         z0_new = SQRT(z0_new**2+rugoro**2)
637         alblw(1:knon) = alb_new(1:knon)         alblw(1:knon) = alb_new(1:knon)
638    
639      else if (nisurf == is_lic) then      else if (nisurf == is_lic) then
640    
641         if (check) write(*,*)'glacier, nisurf = ',nisurf         if (check) write(*, *)'glacier, nisurf = ', nisurf
642    
643         if (.not. allocated(run_off_lic)) then         if (.not. allocated(run_off_lic)) then
644            allocate(run_off_lic(knon), stat = error)            allocate(run_off_lic(knon), stat = error)
645            if (error /= 0) then            if (error /= 0) then
646               abort_message='Pb allocation run_off_lic'               abort_message='Pb allocation run_off_lic'
647               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
648            endif            endif
649            run_off_lic = 0.            run_off_lic = 0.
650         endif         endif
651         !  
652         ! Surface "glacier continentaux" appel a l'interface avec le sol         ! Surface "glacier continentaux" appel a l'interface avec le sol
653         !  
654         IF (soil_model) THEN         IF (soil_model) THEN
655            CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil,soilcap, soilflux)            CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux)
656            cal(1:knon) = RCPD / soilcap(1:knon)            cal(1:knon) = RCPD / soilcap(1:knon)
657            radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)            radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
658         ELSE         ELSE
# Line 701  CONTAINS Line 675  CONTAINS
675              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &              &   radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
676              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &              &   petAcoef, peqAcoef, petBcoef, peqBcoef, &
677              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &              &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
678              &   fqcalving,ffonte, run_off_lic_0)              &   fqcalving, ffonte, run_off_lic_0)
679    
680         ! passage du run-off des glaciers calcule dans fonte_neige au coupleur         ! passage du run-off des glaciers calcule dans fonte_neige au coupleur
681         bidule=0.         bidule=0.
682         bidule(1:knon)= run_off_lic(1:knon)             bidule(1:knon)= run_off_lic(1:knon)    
683         call gath2cpl(bidule, tmp_rlic, klon, knon,iim,jjm,knindex)         call gath2cpl(bidule, tmp_rlic, klon, knon, iim, jjm, knindex)
684         !  
685         ! calcul albedo         ! calcul albedo
686         !  
687         CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))           CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)  
688         WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.         WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
689         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)))
690         alb_new(1 : knon)  = alb_neig(1 : knon)*zfra(1:knon) + &         alb_new(1 : knon)  = alb_neig(1 : knon)*zfra(1:knon) + &
691              &                     0.6 * (1.0-zfra(1:knon))              &                     0.6 * (1.0-zfra(1:knon))
692         !  
693         !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"         !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
694         !       alb_new(1 : knon)  = 0.6 !IM cf FH/GK         !       alb_new(1 : knon)  = 0.6 !IM cf FH/GK
695         !       alb_new(1 : knon)  = 0.82         !       alb_new(1 : knon)  = 0.82
# Line 724  CONTAINS Line 698  CONTAINS
698         !IM: KstaTER0.77 & LMD_ARMIP6             !IM: KstaTER0.77 & LMD_ARMIP6    
699         alb_new(1 : knon)  = 0.77         alb_new(1 : knon)  = 0.77
700    
701         !  
702         ! Rugosite         ! Rugosite
703         !  
704         z0_new = rugoro         z0_new = rugoro
705         !  
706         ! Remplissage des pourcentages de surface         ! Remplissage des pourcentages de surface
707         !  
708         pctsrf_new(:,nisurf) = pctsrf(:,nisurf)         pctsrf_new(:, nisurf) = pctsrf(:, nisurf)
709    
710         alblw(1:knon) = alb_new(1:knon)         alblw(1:knon) = alb_new(1:knon)
711      else      else
712         write(*,*)'Index surface = ',nisurf         write(*, *)'Index surface = ', nisurf
713         abort_message = 'Index surface non valable'         abort_message = 'Index surface non valable'
714         call abort_gcm(modname,abort_message,1)         call abort_gcm(modname, abort_message, 1)
715      endif      endif
716    
717    END SUBROUTINE interfsurf_hq    END SUBROUTINE interfsurf_hq
# Line 747  CONTAINS Line 721  CONTAINS
721    SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &    SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &
722         & radsol, fluxo, fluxg, pctsrf, &         & radsol, fluxo, fluxg, pctsrf, &
723         & tslab, seaice, pctsrf_slab)         & tslab, seaice, pctsrf_slab)
724      !  
725      ! 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
726      ! 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
727      ! la glace de mer pour un "slab" ocean de 50m      ! la glace de mer pour un "slab" ocean de 50m
728      !  
729      ! I. Musat 04.02.2005      ! I. Musat 04.02.2005
730      !  
731      ! input:      ! input:
732      !   klon         nombre total de points de grille      !   klon         nombre total de points de grille
733      !   debut        logical: 1er appel a la physique      !   debut        logical: 1er appel a la physique
# Line 768  CONTAINS Line 742  CONTAINS
742      !   tslab        temperature de l'ocean libre      !   tslab        temperature de l'ocean libre
743      !   seaice       glace de mer (kg/m2)      !   seaice       glace de mer (kg/m2)
744      !   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab      !   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
745      !  
746      use indicesol      use indicesol
747      use clesphys      use clesphys
748      use abort_gcm_m, only: abort_gcm      use abort_gcm_m, only: abort_gcm
# Line 788  CONTAINS Line 762  CONTAINS
762      real, dimension(klon), intent(INOUT) :: tslab      real, dimension(klon), intent(INOUT) :: tslab
763      real, dimension(klon), intent(INOUT)        :: seaice ! glace de mer (kg/m2)      real, dimension(klon), intent(INOUT)        :: seaice ! glace de mer (kg/m2)
764      real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_slab      real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_slab
765      !  
766      ! Variables locales :      ! Variables locales :
767      INTEGER, save :: lmt_pas, julien, idayvrai      INTEGER, save :: lmt_pas, julien, idayvrai
768      REAL, parameter :: unjour=86400.      REAL, parameter :: unjour=86400.
769      real, allocatable, dimension(:), save :: tmp_tslab, tmp_seaice      real, allocatable, dimension(:), save :: tmp_tslab, tmp_seaice
770      REAL, allocatable, dimension(:), save :: slab_bils      REAL, allocatable, dimension(:), save :: slab_bils
771      REAL, allocatable, dimension(:), save :: lmt_bils      REAL, allocatable, dimension(:), save :: lmt_bils
772      logical,save              :: check = .false.      logical, save              :: check = .false.
773      !  
774      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)
775      REAL, parameter :: cbing=0.334e+05        ! J/kg      REAL, parameter :: cbing=0.334e+05        ! J/kg
776      real, dimension(klon)                 :: siceh !hauteur de la glace de mer (m)      real, dimension(klon)                 :: siceh !hauteur de la glace de mer (m)
777      INTEGER :: i      INTEGER :: i
778      integer :: sum_error, error      integer :: sum_error, error
779      REAL :: zz, za, zb      REAL :: zz, za, zb
780      !  
781      character (len = 80) :: abort_message      character (len = 80) :: abort_message
782      character (len = 20) :: modname = 'interfoce_slab'      character (len = 20) :: modname = 'interfoce_slab'
783      !  
784      julien = MOD(ijour,360)      julien = MOD(ijour, 360)
785      sum_error = 0      sum_error = 0
786      IF (debut) THEN      IF (debut) THEN
787         allocate(slab_bils(klon), stat = error); sum_error = sum_error + error         allocate(slab_bils(klon), stat = error); sum_error = sum_error + error
# Line 815  CONTAINS Line 789  CONTAINS
789         allocate(tmp_tslab(klon), stat = error); sum_error = sum_error + error         allocate(tmp_tslab(klon), stat = error); sum_error = sum_error + error
790         allocate(tmp_seaice(klon), stat = error); sum_error = sum_error + error         allocate(tmp_seaice(klon), stat = error); sum_error = sum_error + error
791         if (sum_error /= 0) then         if (sum_error /= 0) then
792            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'
793            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
794         endif         endif
795         tmp_tslab=tslab         tmp_tslab=tslab
796         tmp_seaice=seaice         tmp_seaice=seaice
797         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
798         !  
799         IF (check) THEN         IF (check) THEN
800            PRINT*,'interfoce_slab klon, debut, itap, dtime, ijour, &            PRINT*, 'interfoce_slab klon, debut, itap, dtime, ijour, &
801                 &          lmt_pas ', klon, debut, itap, dtime, ijour, &                 &          lmt_pas ', klon, debut, itap, dtime, ijour, &
802                 &          lmt_pas                 &          lmt_pas
803         ENDIF !check         ENDIF !check
804         !  
805         PRINT*, '************************'         PRINT*, '************************'
806         PRINT*, 'SLAB OCEAN est actif, prenez precautions !'         PRINT*, 'SLAB OCEAN est actif, prenez precautions !'
807         PRINT*, '************************'         PRINT*, '************************'
808         !  
809         ! a mettre un slab_bils aussi en force !!!         ! a mettre un slab_bils aussi en force !!!
810         !  
811         DO i = 1, klon         DO i = 1, klon
812            slab_bils(i) = 0.0            slab_bils(i) = 0.0
813         ENDDO         ENDDO
814         !  
815      ENDIF !debut      ENDIF !debut
816      pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)      pctsrf_slab(1:klon, 1:nbsrf) = pctsrf(1:klon, 1:nbsrf)
817      !  
818      ! 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
819      !  
820      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
821         idayvrai = ijour         idayvrai = ijour
822         CALL condsurf(julien,idayvrai, lmt_bils)         CALL condsurf(julien, idayvrai, lmt_bils)
823      ENDIF !(MOD(itap-1,lmt_pas) .EQ. 0) THEN      ENDIF !(MOD(itap-1, lmt_pas) .EQ. 0) THEN
824    
825      DO i = 1, klon      DO i = 1, klon
826         IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &         IF((pctsrf_slab(i, is_oce).GT.epsfra).OR. &
827              &  (pctsrf_slab(i,is_sic).GT.epsfra)) THEN              &  (pctsrf_slab(i, is_sic).GT.epsfra)) THEN
828            !  
829            ! fabriquer de la glace si congelation atteinte:            ! fabriquer de la glace si congelation atteinte:
830            !  
831            IF (tmp_tslab(i).LT.(RTT-1.8)) THEN            IF (tmp_tslab(i).LT.(RTT-1.8)) THEN
832               zz =  (RTT-1.8)-tmp_tslab(i)               zz =  (RTT-1.8)-tmp_tslab(i)
833               tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz               tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz
834               seaice(i) = tmp_seaice(i)               seaice(i) = tmp_seaice(i)
835               tmp_tslab(i) = RTT-1.8               tmp_tslab(i) = RTT-1.8
836            ENDIF            ENDIF
837            !  
838            ! faire fondre de la glace si temperature est superieure a 0:            ! faire fondre de la glace si temperature est superieure a 0:
839            !  
840            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
841               zz = cyang/cbing * (tmp_tslab(i)-RTT)               zz = cyang/cbing * (tmp_tslab(i)-RTT)
842               zz = MIN(zz,tmp_seaice(i))               zz = MIN(zz, tmp_seaice(i))
843               tmp_seaice(i) = tmp_seaice(i) - zz               tmp_seaice(i) = tmp_seaice(i) - zz
844               seaice(i) = tmp_seaice(i)               seaice(i) = tmp_seaice(i)
845               tmp_tslab(i) = tmp_tslab(i) - zz*cbing/cyang               tmp_tslab(i) = tmp_tslab(i) - zz*cbing/cyang
846            ENDIF            ENDIF
847            !  
848            ! limiter la glace de mer a 10 metres (10000 kg/m2)            ! limiter la glace de mer a 10 metres (10000 kg/m2)
849            !  
850            IF(tmp_seaice(i).GT.45.) THEN            IF(tmp_seaice(i).GT.45.) THEN
851               tmp_seaice(i) = MIN(tmp_seaice(i),10000.0)               tmp_seaice(i) = MIN(tmp_seaice(i), 10000.0)
852            ELSE            ELSE
853               tmp_seaice(i) = 0.               tmp_seaice(i) = 0.
854            ENDIF            ENDIF
855            seaice(i) = tmp_seaice(i)            seaice(i) = tmp_seaice(i)
856            siceh(i)=tmp_seaice(i)/1000. !en metres            siceh(i)=tmp_seaice(i)/1000. !en metres
857            !  
858            ! determiner la nature du sol (glace de mer ou ocean libre):            ! determiner la nature du sol (glace de mer ou ocean libre):
859            !  
860            ! on fait dependre la fraction de seaice "pctsrf(i,is_sic)"            ! on fait dependre la fraction de seaice "pctsrf(i, is_sic)"
861            ! de l'epaisseur de seaice :            ! de l'epaisseur de seaice :
862            ! 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
863            ! 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
864            !  
865            pctsrf_slab(i,is_sic)=MIN(siceh(i)/0.20, &            pctsrf_slab(i, is_sic)=MIN(siceh(i)/0.20, &
866                 &                      1.-(pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)))                 &                      1.-(pctsrf_slab(i, is_ter)+pctsrf_slab(i, is_lic)))
867            pctsrf_slab(i,is_oce)=1.0 - &            pctsrf_slab(i, is_oce)=1.0 - &
868                 &      (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))
869         ENDIF !pctsrf         ENDIF !pctsrf
870      ENDDO      ENDDO
871      !  
872      ! Calculer le bilan du flux de chaleur au sol :      ! Calculer le bilan du flux de chaleur au sol :
873      !  
874      DO i = 1, klon      DO i = 1, klon
875         za = radsol(i) + fluxo(i)         za = radsol(i) + fluxo(i)
876         zb = fluxg(i)         zb = fluxg(i)
877         IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &         IF((pctsrf_slab(i, is_oce).GT.epsfra).OR. &
878              &   (pctsrf_slab(i,is_sic).GT.epsfra)) THEN              &   (pctsrf_slab(i, is_sic).GT.epsfra)) THEN
879            slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i,is_oce) &            slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i, is_oce) &
880                 &             +zb*pctsrf_slab(i,is_sic))/ FLOAT(lmt_pas)                 &             +zb*pctsrf_slab(i, is_sic))/ FLOAT(lmt_pas)
881         ENDIF         ENDIF
882      ENDDO !klon      ENDDO !klon
883      !  
884      ! calcul tslab      ! calcul tslab
885      !  
886      IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee      IF (MOD(itap, lmt_pas).EQ.0) THEN !fin de journee
887         DO i = 1, klon         DO i = 1, klon
888            IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &            IF ((pctsrf_slab(i, is_oce).GT.epsfra).OR. &
889                 &    (pctsrf_slab(i,is_sic).GT.epsfra)) THEN                 &    (pctsrf_slab(i, is_sic).GT.epsfra)) THEN
890               tmp_tslab(i) = tmp_tslab(i) + &               tmp_tslab(i) = tmp_tslab(i) + &
891                    & (slab_bils(i)-lmt_bils(i)) &                    & (slab_bils(i)-lmt_bils(i)) &
892                    &                         /cyang*unjour                    &                         /cyang*unjour
# Line 920  CONTAINS Line 894  CONTAINS
894               slab_bils(i) = 0.               slab_bils(i) = 0.
895            ENDIF !pctsrf            ENDIF !pctsrf
896         ENDDO !klon         ENDDO !klon
897      ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN      ENDIF !(MOD(itap, lmt_pas).EQ.0) THEN
898      !  
899      tslab = tmp_tslab      tslab = tmp_tslab
900      seaice = tmp_seaice      seaice = tmp_seaice
901    END SUBROUTINE interfoce_slab    END SUBROUTINE interfoce_slab
# Line 933  CONTAINS Line 907  CONTAINS
907         & debut,  &         & debut,  &
908         & lmt_sst, pctsrf_new)         & lmt_sst, pctsrf_new)
909    
910      ! Cette routine sert d'interface entre le modele atmospherique et un fichier      ! Cette routine sert d'interface entre le modele atmospherique et
911      ! de conditions aux limites      ! un fichier de conditions aux limites
912      !  
913      ! 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  
     !  
914    
915      use abort_gcm_m, only: abort_gcm      use abort_gcm_m, only: abort_gcm
916      use indicesol      use indicesol
917    
918      ! Parametres d'entree      integer, intent(IN) :: itime ! numero du pas de temps courant
919      integer, intent(IN) :: itime      real   , intent(IN) :: dtime ! pas de temps de la physique (en s)
920      real   , intent(IN) :: dtime      integer, intent(IN) :: jour ! jour a lire dans l'annee
921      integer, intent(IN) :: jour      integer, intent(IN) :: nisurf ! index de la surface a traiter (1 = sol continental)
922      integer, intent(IN) :: nisurf      integer, intent(IN) :: knon ! nombre de points dans le domaine a traiter
923      integer, intent(IN) :: knon      integer, intent(IN) :: klon ! taille de la grille
924      integer, intent(IN) :: klon      integer, dimension(klon), intent(in) :: knindex ! index des points de la surface a traiter
925      integer, dimension(klon), intent(in) :: knindex      logical, intent(IN) :: debut ! logical: 1er appel a la physique (initialisation)
     logical, intent(IN) :: debut  
926    
927      ! Parametres de sortie      ! Parametres de sortie
928        ! output:
929        !   lmt_sst      SST lues dans le fichier de CL
930        !   pctsrf_new   sous-maille fractionnelle
931      real, intent(out), dimension(klon) :: lmt_sst      real, intent(out), dimension(klon) :: lmt_sst
932      real, intent(out), dimension(klon,nbsrf) :: pctsrf_new      real, intent(out), dimension(klon, nbsrf) :: pctsrf_new
933    
934      ! Variables locales      ! Variables locales
935      integer     :: ii      integer     :: ii
936      INTEGER,save :: lmt_pas     ! frequence de lecture des conditions limites      INTEGER, save :: lmt_pas     ! frequence de lecture des conditions limites
937      ! (en pas de physique)      ! (en pas de physique)
938      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
939      ! lu pour une surface precedente      ! lu pour une surface precedente
940      integer,save :: jour_lu      integer, save :: jour_lu
941      integer      :: ierr      integer      :: ierr
942      character (len = 20) :: modname = 'interfoce_lim'      character (len = 20) :: modname = 'interfoce_lim'
943      character (len = 80) :: abort_message      character (len = 80) :: abort_message
# Line 984  CONTAINS Line 945  CONTAINS
945      logical, save     :: check = .FALSE.      logical, save     :: check = .FALSE.
946      ! Champs lus dans le fichier de CL      ! Champs lus dans le fichier de CL
947      real, allocatable , save, dimension(:) :: sst_lu, rug_lu, nat_lu      real, allocatable , save, dimension(:) :: sst_lu, rug_lu, nat_lu
948      real, allocatable , save, dimension(:,:) :: pct_tmp      real, allocatable , save, dimension(:, :) :: pct_tmp
949      !  
950      ! quelques variables pour netcdf      ! quelques variables pour netcdf
951      !  
952      include "netcdf.inc"      include "netcdf.inc"
953      integer              :: nid, nvarid      integer              :: nid, nvarid
954      integer, dimension(2) :: start, epais      integer, dimension(2) :: start, epais
955      !  
956      ! Fin déclaration      ! --------------------------------------------------
     !  
957    
958      if (debut .and. .not. allocated(sst_lu)) then      if (debut .and. .not. allocated(sst_lu)) then
959         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
960         jour_lu = jour - 1         jour_lu = jour - 1
961         allocate(sst_lu(klon))         allocate(sst_lu(klon))
962         allocate(nat_lu(klon))         allocate(nat_lu(klon))
963         allocate(pct_tmp(klon,nbsrf))         allocate(pct_tmp(klon, nbsrf))
964      endif      endif
965    
966      if ((jour - jour_lu) /= 0) deja_lu = .false.      if ((jour - jour_lu) /= 0) deja_lu = .false.
967    
968      if (check) write(*,*)modname,' :: jour, jour_lu, deja_lu', jour, jour_lu, &      if (check) write(*, *)modname, ' :: jour, jour_lu, deja_lu', jour, jour_lu, &
969           deja_lu           deja_lu
970      if (check) write(*,*)modname,' :: itime, lmt_pas ', itime, lmt_pas,dtime      if (check) write(*, *)modname, ' :: itime, lmt_pas ', itime, lmt_pas, dtime
971    
972      ! Tester d'abord si c'est le moment de lire le fichier      ! Tester d'abord si c'est le moment de lire le fichier
973      if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then      if (mod(itime-1, lmt_pas) == 0 .and. .not. deja_lu) then
974         !  
975         ! Ouverture du fichier         ! Ouverture du fichier
976         !  
977         ierr = NF_OPEN ('limit.nc', NF_NOWRITE,nid)         ierr = NF_OPEN ('limit.nc', NF_NOWRITE, nid)
978         if (ierr.NE.NF_NOERR) then         if (ierr.NE.NF_NOERR) then
979            abort_message &            abort_message &
980                 = 'Pb d''ouverture du fichier de conditions aux limites'                 = 'Pb d''ouverture du fichier de conditions aux limites'
981            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
982         endif         endif
983         !  
984         ! La tranche de donnees a lire:         ! La tranche de donnees a lire:
985         !  
986         start(1) = 1         start(1) = 1
987         start(2) = jour         start(2) = jour
988         epais(1) = klon         epais(1) = klon
989         epais(2) = 1         epais(2) = 1
990         !  
991         if (newlmt) then         if (newlmt) then
992            !  
993            ! Fraction "ocean"            ! Fraction "ocean"
994            !  
995            ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)            ierr = NF_INQ_VARID(nid, 'FOCE', nvarid)
996            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
997               abort_message = 'Le champ <FOCE> est absent'               abort_message = 'Le champ <FOCE> 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_oce))            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_oce))
1001            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1002               abort_message = 'Lecture echouee pour <FOCE>'               abort_message = 'Lecture echouee pour <FOCE>'
1003               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1004            endif            endif
1005            !  
1006            ! Fraction "glace de mer"            ! Fraction "glace de mer"
1007            !  
1008            ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)            ierr = NF_INQ_VARID(nid, 'FSIC', nvarid)
1009            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1010               abort_message = 'Le champ <FSIC> est absent'               abort_message = 'Le champ <FSIC> 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_sic))            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_sic))
1014            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1015               abort_message = 'Lecture echouee pour <FSIC>'               abort_message = 'Lecture echouee pour <FSIC>'
1016               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1017            endif            endif
1018            !  
1019            ! Fraction "terre"            ! Fraction "terre"
1020            !  
1021            ierr = NF_INQ_VARID(nid, 'FTER', nvarid)            ierr = NF_INQ_VARID(nid, 'FTER', nvarid)
1022            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1023               abort_message = 'Le champ <FTER> est absent'               abort_message = 'Le champ <FTER> 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_ter))            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_ter))
1027            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1028               abort_message = 'Lecture echouee pour <FTER>'               abort_message = 'Lecture echouee pour <FTER>'
1029               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1030            endif            endif
1031            !  
1032            ! Fraction "glacier terre"            ! Fraction "glacier terre"
1033            !  
1034            ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)            ierr = NF_INQ_VARID(nid, 'FLIC', nvarid)
1035            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1036               abort_message = 'Le champ <FLIC> est absent'               abort_message = 'Le champ <FLIC> 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,pct_tmp(1,is_lic))            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, pct_tmp(1, is_lic))
1040            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1041               abort_message = 'Lecture echouee pour <FLIC>'               abort_message = 'Lecture echouee pour <FLIC>'
1042               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1043            endif            endif
1044            !  
1045         else  ! on en est toujours a rnatur         else  ! on en est toujours a rnatur
1046            !  
1047            ierr = NF_INQ_VARID(nid, 'NAT', nvarid)            ierr = NF_INQ_VARID(nid, 'NAT', nvarid)
1048            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1049               abort_message = 'Le champ <NAT> est absent'               abort_message = 'Le champ <NAT> est absent'
1050               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1051            endif            endif
1052            ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, nat_lu)            ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, nat_lu)
1053            if (ierr /= NF_NOERR) then            if (ierr /= NF_NOERR) then
1054               abort_message = 'Lecture echouee pour <NAT>'               abort_message = 'Lecture echouee pour <NAT>'
1055               call abort_gcm(modname,abort_message,1)               call abort_gcm(modname, abort_message, 1)
1056            endif            endif
1057            !  
1058            ! Remplissage des fractions de surface            ! Remplissage des fractions de surface
1059            ! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice            ! nat = 0, 1, 2, 3 pour ocean, terre, glacier, seaice
1060            !  
1061            pct_tmp = 0.0            pct_tmp = 0.0
1062            do ii = 1, klon            do ii = 1, klon
1063               pct_tmp(ii,nint(nat_lu(ii)) + 1) = 1.               pct_tmp(ii, nint(nat_lu(ii)) + 1) = 1.
1064            enddo            enddo
1065    
1066            !  
1067            !  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
1068            !  
1069            pctsrf_new = pct_tmp            pctsrf_new = pct_tmp
1070            pctsrf_new (:,2)= pct_tmp (:,1)            pctsrf_new (:, 2)= pct_tmp (:, 1)
1071            pctsrf_new (:,1)= pct_tmp (:,2)            pctsrf_new (:, 1)= pct_tmp (:, 2)
1072            pct_tmp = pctsrf_new            pct_tmp = pctsrf_new
1073         endif ! fin test sur newlmt         endif ! fin test sur newlmt
1074         !  
1075         ! Lecture SST         ! Lecture SST
1076         !  
1077         ierr = NF_INQ_VARID(nid, 'SST', nvarid)         ierr = NF_INQ_VARID(nid, 'SST', nvarid)
1078         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1079            abort_message = 'Le champ <SST> est absent'            abort_message = 'Le champ <SST> est absent'
1080            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1081         endif         endif
1082         ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, sst_lu)         ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, sst_lu)
1083         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1084            abort_message = 'Lecture echouee pour <SST>'            abort_message = 'Lecture echouee pour <SST>'
1085            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1086         endif         endif
1087    
1088         !  
1089         ! Fin de lecture         ! Fin de lecture
1090         !  
1091         ierr = NF_CLOSE(nid)         ierr = NF_CLOSE(nid)
1092         deja_lu = .true.         deja_lu = .true.
1093         jour_lu = jour         jour_lu = jour
1094      endif      endif
1095      !  
1096      ! Recopie des variables dans les champs de sortie      ! Recopie des variables dans les champs de sortie
1097      !  
1098      lmt_sst = 999999999.      lmt_sst = 999999999.
1099      do ii = 1, knon      do ii = 1, knon
1100         lmt_sst(ii) = sst_lu(knindex(ii))         lmt_sst(ii) = sst_lu(knindex(ii))
1101      enddo      enddo
1102    
1103      pctsrf_new(:,is_oce) = pct_tmp(:,is_oce)      pctsrf_new(:, is_oce) = pct_tmp(:, is_oce)
1104      pctsrf_new(:,is_sic) = pct_tmp(:,is_sic)      pctsrf_new(:, is_sic) = pct_tmp(:, is_sic)
1105    
1106    END SUBROUTINE interfoce_lim    END SUBROUTINE interfoce_lim
1107    
# Line 1154  CONTAINS Line 1114  CONTAINS
1114    
1115      ! Cette routine sert d'interface entre le modèle atmosphérique et      ! Cette routine sert d'interface entre le modèle atmosphérique et
1116      ! un fichier de conditions aux limites.      ! un fichier de conditions aux limites.
1117      !  
1118      ! L. Fairhead 02/2000      ! L. Fairhead 02/2000
1119      !  
1120        use abort_gcm_m, only: abort_gcm
1121    
1122        ! Parametres d'entree
1123      ! input:      ! input:
1124      !   itime        numero du pas de temps courant      !   itime        numero du pas de temps courant
1125      !   dtime        pas de temps de la physique (en s)      !   dtime        pas de temps de la physique (en s)
# Line 1166  CONTAINS Line 1129  CONTAINS
1129      !   knindex      index des points de la surface a traiter      !   knindex      index des points de la surface a traiter
1130      !   klon         taille de la grille      !   klon         taille de la grille
1131      !   debut        logical: 1er appel a la physique (initialisation)      !   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  
     !  
   
     use abort_gcm_m, only: abort_gcm  
   
     ! Parametres d'entree  
1132      integer, intent(IN) :: itime      integer, intent(IN) :: itime
1133      real   , intent(IN) :: dtime      real   , intent(IN) :: dtime
1134      integer, intent(IN) :: jour      integer, intent(IN) :: jour
# Line 1187  CONTAINS Line 1139  CONTAINS
1139      logical, intent(IN) :: debut      logical, intent(IN) :: debut
1140    
1141      ! Parametres de sortie      ! Parametres de sortie
1142        ! output:
1143        !   lmt_sst      SST lues dans le fichier de CL
1144        !   lmt_alb      Albedo lu
1145        !   lmt_rug      longueur de rugosité lue
1146        !   pctsrf_new   sous-maille fractionnelle
1147      real, intent(out), dimension(klon) :: lmt_alb      real, intent(out), dimension(klon) :: lmt_alb
1148      real, intent(out), dimension(klon) :: lmt_rug      real, intent(out), dimension(klon) :: lmt_rug
1149    
1150      ! Variables locales      ! Variables locales
1151      integer     :: ii      integer     :: ii
1152      integer,save :: lmt_pas     ! frequence de lecture des conditions limites      integer, save :: lmt_pas     ! frequence de lecture des conditions limites
1153      ! (en pas de physique)      ! (en pas de physique)
1154      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
1155      ! lu pour une surface precedente      ! lu pour une surface precedente
1156      integer,save :: jour_lu_sur      integer, save :: jour_lu_sur
1157      integer      :: ierr      integer      :: ierr
1158      character (len = 20) :: modname = 'interfsur_lim'      character (len = 20) :: modname = 'interfsur_lim'
1159      character (len = 80) :: abort_message      character (len = 80) :: abort_message
1160      logical,save     :: newlmt = .false.      logical, save     :: newlmt = .false.
1161      logical,save     :: check = .false.      logical, save     :: check = .false.
1162      ! Champs lus dans le fichier de CL      ! Champs lus dans le fichier de CL
1163      real, allocatable , save, dimension(:) :: alb_lu, rug_lu      real, allocatable , save, dimension(:) :: alb_lu, rug_lu
1164      !  
1165      ! quelques variables pour netcdf      ! quelques variables pour netcdf
1166      !  
1167      include "netcdf.inc"      include "netcdf.inc"
1168      integer ,save             :: nid, nvarid      integer , save             :: nid, nvarid
1169      integer, dimension(2),save :: start, epais      integer, dimension(2), save :: start, epais
1170      !  
1171      ! Fin déclaration      !------------------------------------------------------------
     !  
1172    
1173      if (debut) then      if (debut) then
1174         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 1223  CONTAINS Line 1179  CONTAINS
1179    
1180      if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false.      if ((jour - jour_lu_sur) /= 0) deja_lu_sur = .false.
1181    
1182      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, &
1183           deja_lu_sur           deja_lu_sur
1184      if (check) write(*,*)modname,':: itime, lmt_pas', itime, lmt_pas      if (check) write(*, *)modname, ':: itime, lmt_pas', itime, lmt_pas
1185    
1186      ! Tester d'abord si c'est le moment de lire le fichier      ! Tester d'abord si c'est le moment de lire le fichier
1187      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
1188         !    
1189         ! Ouverture du fichier         ! Ouverture du fichier
1190         !    
1191         ierr = NF_OPEN ('limit.nc', NF_NOWRITE,nid)         ierr = NF_OPEN ('limit.nc', NF_NOWRITE, nid)
1192         if (ierr.NE.NF_NOERR) then         if (ierr.NE.NF_NOERR) then
1193            abort_message &            abort_message &
1194                 = 'Pb d''ouverture du fichier de conditions aux limites'                 = 'Pb d''ouverture du fichier de conditions aux limites'
1195            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1196         endif         endif
1197         !  
1198         ! La tranche de donnees a lire:         ! La tranche de donnees a lire:
1199    
1200         start(1) = 1         start(1) = 1
1201         start(2) = jour         start(2) = jour
1202         epais(1) = klon         epais(1) = klon
1203         epais(2) = 1         epais(2) = 1
1204         !    
1205         ! Lecture Albedo         ! Lecture Albedo
1206         !    
1207         ierr = NF_INQ_VARID(nid, 'ALB', nvarid)         ierr = NF_INQ_VARID(nid, 'ALB', nvarid)
1208         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1209            abort_message = 'Le champ <ALB> est absent'            abort_message = 'Le champ <ALB> 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, alb_lu)         ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, alb_lu)
1213         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1214            abort_message = 'Lecture echouee pour <ALB>'            abort_message = 'Lecture echouee pour <ALB>'
1215            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1216         endif         endif
1217         !    
1218         ! Lecture rugosité         ! Lecture rugosité
1219         !    
1220         ierr = NF_INQ_VARID(nid, 'RUG', nvarid)         ierr = NF_INQ_VARID(nid, 'RUG', nvarid)
1221         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1222            abort_message = 'Le champ <RUG> est absent'            abort_message = 'Le champ <RUG> est absent'
1223            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1224         endif         endif
1225         ierr = NF_GET_VARA_REAL(nid,nvarid,start,epais, rug_lu)         ierr = NF_GET_VARA_REAL(nid, nvarid, start, epais, rug_lu)
1226         if (ierr /= NF_NOERR) then         if (ierr /= NF_NOERR) then
1227            abort_message = 'Lecture echouee pour <RUG>'            abort_message = 'Lecture echouee pour <RUG>'
1228            call abort_gcm(modname,abort_message,1)            call abort_gcm(modname, abort_message, 1)
1229         endif         endif
1230    
1231         !    
1232         ! Fin de lecture         ! Fin de lecture
1233         !    
1234         ierr = NF_CLOSE(nid)         ierr = NF_CLOSE(nid)
1235         deja_lu_sur = .true.         deja_lu_sur = .true.
1236         jour_lu_sur = jour         jour_lu_sur = jour
1237      endif      endif
1238      !  
1239      ! Recopie des variables dans les champs de sortie      ! Recopie des variables dans les champs de sortie
1240      !  
1241  !!$  lmt_alb(:) = 0.0  !!$  lmt_alb = 0.0
1242  !!$  lmt_rug(:) = 0.0  !!$  lmt_rug = 0.0
1243      lmt_alb(:) = 999999.      lmt_alb = 999999.
1244      lmt_rug(:) = 999999.      lmt_rug = 999999.
1245      DO ii = 1, knon      DO ii = 1, knon
1246         lmt_alb(ii) = alb_lu(knindex(ii))         lmt_alb(ii) = alb_lu(knindex(ii))
1247         lmt_rug(ii) = rug_lu(knindex(ii))         lmt_rug(ii) = rug_lu(knindex(ii))
# Line 1304  CONTAINS Line 1260  CONTAINS
1260    
1261      ! 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
1262      ! une temperature de surface (au cas ou ok_veget = false)      ! une temperature de surface (au cas ou ok_veget = false)
1263      !  
1264      ! L. Fairhead 4/2000      ! L. Fairhead 4/2000
1265      !  
1266      ! input:      ! input:
1267      !   knon         nombre de points a traiter      !   knon         nombre de points a traiter
1268      !   nisurf       surface a traiter      !   nisurf       surface a traiter
# Line 1326  CONTAINS Line 1282  CONTAINS
1282      !   peqBcoef     coeff. B de la resolution de la CL pour q      !   peqBcoef     coeff. B de la resolution de la CL pour q
1283      !   radsol       rayonnement net aus sol (LW + SW)      !   radsol       rayonnement net aus sol (LW + SW)
1284      !   dif_grnd     coeff. diffusion vers le sol profond      !   dif_grnd     coeff. diffusion vers le sol profond
1285      !  
1286      ! output:      ! output:
1287      !   tsurf_new    temperature au sol      !   tsurf_new    temperature au sol
1288      !   qsurf        humidite de l'air au dessus du sol      !   qsurf        humidite de l'air au dessus du sol
# Line 1334  CONTAINS Line 1290  CONTAINS
1290      !   fluxlat      flux de chaleur latente      !   fluxlat      flux de chaleur latente
1291      !   dflux_s      derivee du flux de chaleur sensible / Ts      !   dflux_s      derivee du flux de chaleur sensible / Ts
1292      !   dflux_l      derivee du flux de chaleur latente  / Ts      !   dflux_l      derivee du flux de chaleur latente  / Ts
1293      !  
1294    
1295      use indicesol      use indicesol
1296      use abort_gcm_m, only: abort_gcm      use abort_gcm_m, only: abort_gcm
# Line 1372  CONTAINS Line 1328  CONTAINS
1328      real, parameter :: t_grnd = 271.35, t_coup = 273.15      real, parameter :: t_grnd = 271.35, t_coup = 273.15
1329      !! PB temporaire en attendant mieux pour le modele de neige      !! PB temporaire en attendant mieux pour le modele de neige
1330      REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)      REAL, parameter :: chasno = 3.334E+05/(2.3867E+06*0.15)
1331      !  
1332      logical, save         :: check = .false.      logical, save         :: check = .false.
1333      character (len = 20)  :: modname = 'calcul_fluxs'      character (len = 20)  :: modname = 'calcul_fluxs'
1334      logical, save         :: fonte_neige = .false.      logical, save         :: fonte_neige = .false.
1335      real, save            :: max_eau_sol = 150.0      real, save            :: max_eau_sol = 150.0
1336      character (len = 80) :: abort_message      character (len = 80) :: abort_message
1337      logical,save         :: first = .true.,second=.false.      logical, save         :: first = .true., second=.false.
1338    
1339      if (check) write(*,*)'Entree ', modname,' surface = ',nisurf      if (check) write(*, *)'Entree ', modname, ' surface = ', nisurf
1340    
1341      IF (check) THEN      IF (check) THEN
1342         WRITE(*,*)' radsol (min, max)' &         WRITE(*, *)' radsol (min, max)' &
1343              &     , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))              &     , MINVAL(radsol(1:knon)), MAXVAL(radsol(1:knon))
1344         !!CALL flush(6)         !!CALL flush(6)
1345      ENDIF      ENDIF
1346    
1347      if (size(coastalflow) /= knon .AND. nisurf == is_ter) then      if (size(coastalflow) /= knon .AND. nisurf == is_ter) then
1348         write(*,*)'Bizarre, le nombre de points continentaux'         write(*, *)'Bizarre, le nombre de points continentaux'
1349         write(*,*)'a change entre deux appels. J''arrete ...'         write(*, *)'a change entre deux appels. J''arrete ...'
1350         abort_message='Pb run_off'         abort_message='Pb run_off'
1351         call abort_gcm(modname,abort_message,1)         call abort_gcm(modname, abort_message, 1)
1352      endif      endif
1353      !  
1354      ! Traitement neige et humidite du sol      ! Traitement neige et humidite du sol
1355      !  
1356  !!$  WRITE(*,*)'test calcul_flux, surface ', nisurf  !!$  WRITE(*, *)'test calcul_flux, surface ', nisurf
1357      !!PB test      !!PB test
1358  !!$    if (nisurf == is_oce) then  !!$    if (nisurf == is_oce) then
1359  !!$      snow = 0.  !!$      snow = 0.
# Line 1410  CONTAINS Line 1366  CONTAINS
1366  !!$    endif  !!$    endif
1367  !!$    IF (nisurf /= is_ter) qsol = max_eau_sol  !!$    IF (nisurf /= is_ter) qsol = max_eau_sol
1368    
1369      !  
1370      ! Initialisation      ! Initialisation
1371      !  
1372      evap = 0.      evap = 0.
1373      fluxsens=0.      fluxsens=0.
1374      fluxlat=0.      fluxlat=0.
1375      dflux_s = 0.      dflux_s = 0.
1376      dflux_l = 0.      dflux_l = 0.
1377      !  
1378      ! zx_qs = qsat en kg/kg      ! zx_qs = qsat en kg/kg
1379      !  
1380      DO i = 1, knon      DO i = 1, knon
1381         zx_pkh(i) = (ps(i)/ps(i))**RKAPPA         zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
1382         IF (thermcep) THEN         IF (thermcep) THEN
1383            zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))            zdelta=MAX(0., SIGN(1., rtt-tsurf(i)))
1384            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
1385            zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))            zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
1386            zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)            zx_qs= r2es * FOEEW(tsurf(i), zdelta)/ps(i)
1387            zx_qs=MIN(0.5,zx_qs)            zx_qs=MIN(0.5, zx_qs)
1388            zcor=1./(1.-retv*zx_qs)            zcor=1./(1.-retv*zx_qs)
1389            zx_qs=zx_qs*zcor            zx_qs=zx_qs*zcor
1390            zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &            zx_dq_s_dh = FOEDE(tsurf(i), zdelta, zcvm5, zx_qs, zcor) &
1391                 &                 /RLVTT / zx_pkh(i)                 &                 /RLVTT / zx_pkh(i)
1392         ELSE         ELSE
1393            IF (tsurf(i).LT.t_coup) THEN            IF (tsurf(i).LT.t_coup) THEN
1394               zx_qs = qsats(tsurf(i)) / ps(i)               zx_qs = qsats(tsurf(i)) / ps(i)
1395               zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &               zx_dq_s_dh = dqsats(tsurf(i), zx_qs)/RLVTT &
1396                    &                    / zx_pkh(i)                    &                    / zx_pkh(i)
1397            ELSE            ELSE
1398               zx_qs = qsatl(tsurf(i)) / ps(i)               zx_qs = qsatl(tsurf(i)) / ps(i)
1399               zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &               zx_dq_s_dh = dqsatl(tsurf(i), zx_qs)/RLVTT &
1400                    &               / zx_pkh(i)                    &               / zx_pkh(i)
1401            ENDIF            ENDIF
1402         ENDIF         ENDIF
# Line 1453  CONTAINS Line 1409  CONTAINS
1409      ENDDO      ENDDO
1410    
1411      ! === Calcul de la temperature de surface ===      ! === Calcul de la temperature de surface ===
1412      !  
1413      ! zx_sl = chaleur latente d'evaporation ou de sublimation      ! zx_sl = chaleur latente d'evaporation ou de sublimation
1414      !  
1415      do i = 1, knon      do i = 1, knon
1416         zx_sl(i) = RLVTT         zx_sl(i) = RLVTT
1417         if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT         if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
# Line 1485  CONTAINS Line 1441  CONTAINS
1441              &                       zx_nh(i) + zx_sl(i) * zx_nq(i)) &                &                       zx_nh(i) + zx_sl(i) * zx_nq(i)) &  
1442              &                     + dtime * dif_grnd(i))              &                     + dtime * dif_grnd(i))
1443    
1444         !  
1445         ! Y'a-t-il fonte de neige?         ! Y'a-t-il fonte de neige?
1446         !  
1447         !    fonte_neige = (nisurf /= is_oce) .AND. &         !    fonte_neige = (nisurf /= is_oce) .AND. &
1448         !     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &         !     & (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
1449         !     & .AND. (tsurf_new(i) >= RTT)         !     & .AND. (tsurf_new(i) >= RTT)
# Line 1519  CONTAINS Line 1475  CONTAINS
1475         & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &         & radsol, dif_grnd, t1lay, q1lay, u1lay, v1lay, &
1476         & petAcoef, peqAcoef, petBcoef, peqBcoef, &         & petAcoef, peqAcoef, petBcoef, peqBcoef, &
1477         & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &         & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
1478         & fqcalving,ffonte,run_off_lic_0)         & fqcalving, ffonte, run_off_lic_0)
1479    
1480      ! 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
1481      ! de sol simplifié      ! de sol simplifié
1482      !  
1483      ! LF 03/2001      ! LF 03/2001
1484      ! input:      ! input:
1485      !   knon         nombre de points a traiter      !   knon         nombre de points a traiter
# Line 1545  CONTAINS Line 1501  CONTAINS
1501      !   peqBcoef     coeff. B de la resolution de la CL pour q      !   peqBcoef     coeff. B de la resolution de la CL pour q
1502      !   radsol       rayonnement net aus sol (LW + SW)      !   radsol       rayonnement net aus sol (LW + SW)
1503      !   dif_grnd     coeff. diffusion vers le sol profond      !   dif_grnd     coeff. diffusion vers le sol profond
1504      !  
1505      ! output:      ! output:
1506      !   tsurf_new    temperature au sol      !   tsurf_new    temperature au sol
1507      !   fluxsens     flux de chaleur sensible      !   fluxsens     flux de chaleur sensible
# Line 1554  CONTAINS Line 1510  CONTAINS
1510      !   dflux_l      derivee du flux de chaleur latente  / Ts      !   dflux_l      derivee du flux de chaleur latente  / Ts
1511      ! in/out:      ! in/out:
1512      !   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
1513      !  
1514    
1515      use indicesol      use indicesol
1516      use YOMCST      use YOMCST
# Line 1606  CONTAINS Line 1562  CONTAINS
1562      !IM cf JLD/ GKtest      !IM cf JLD/ GKtest
1563      REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)      REAL, parameter :: chaice = 3.334E+05/(2.3867E+06*0.15)
1564      ! fin GKtest      ! fin GKtest
1565      !  
1566      logical, save         :: check = .FALSE.      logical, save         :: check = .FALSE.
1567      character (len = 20)  :: modname = 'fonte_neige'      character (len = 20)  :: modname = 'fonte_neige'
1568      logical, save         :: neige_fond = .false.      logical, save         :: neige_fond = .false.
1569      real, save            :: max_eau_sol = 150.0      real, save            :: max_eau_sol = 150.0
1570      character (len = 80) :: abort_message      character (len = 80) :: abort_message
1571      logical,save         :: first = .true.,second=.false.      logical, save         :: first = .true., second=.false.
1572      real                 :: coeff_rel      real                 :: coeff_rel
1573    
1574      if (check) write(*,*)'Entree ', modname,' surface = ',nisurf      if (check) write(*, *)'Entree ', modname, ' surface = ', nisurf
1575    
1576      ! Initialisations      ! Initialisations
1577      coeff_rel = dtime/(tau_calv * rday)      coeff_rel = dtime/(tau_calv * rday)
1578      bil_eau_s(:) = 0.      bil_eau_s = 0.
1579      DO i = 1, knon      DO i = 1, knon
1580         zx_pkh(i) = (ps(i)/ps(i))**RKAPPA         zx_pkh(i) = (ps(i)/ps(i))**RKAPPA
1581         IF (thermcep) THEN         IF (thermcep) THEN
1582            zdelta=MAX(0.,SIGN(1.,rtt-tsurf(i)))            zdelta=MAX(0., SIGN(1., rtt-tsurf(i)))
1583            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta            zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
1584            zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))            zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*q1lay(i))
1585            zx_qs= r2es * FOEEW(tsurf(i),zdelta)/ps(i)            zx_qs= r2es * FOEEW(tsurf(i), zdelta)/ps(i)
1586            zx_qs=MIN(0.5,zx_qs)            zx_qs=MIN(0.5, zx_qs)
1587            zcor=1./(1.-retv*zx_qs)            zcor=1./(1.-retv*zx_qs)
1588            zx_qs=zx_qs*zcor            zx_qs=zx_qs*zcor
1589            zx_dq_s_dh = FOEDE(tsurf(i),zdelta,zcvm5,zx_qs,zcor) &            zx_dq_s_dh = FOEDE(tsurf(i), zdelta, zcvm5, zx_qs, zcor) &
1590                 &                 /RLVTT / zx_pkh(i)                 &                 /RLVTT / zx_pkh(i)
1591         ELSE         ELSE
1592            IF (tsurf(i).LT.t_coup) THEN            IF (tsurf(i).LT.t_coup) THEN
1593               zx_qs = qsats(tsurf(i)) / ps(i)               zx_qs = qsats(tsurf(i)) / ps(i)
1594               zx_dq_s_dh = dqsats(tsurf(i),zx_qs)/RLVTT &               zx_dq_s_dh = dqsats(tsurf(i), zx_qs)/RLVTT &
1595                    &                    / zx_pkh(i)                    &                    / zx_pkh(i)
1596            ELSE            ELSE
1597               zx_qs = qsatl(tsurf(i)) / ps(i)               zx_qs = qsatl(tsurf(i)) / ps(i)
1598               zx_dq_s_dh = dqsatl(tsurf(i),zx_qs)/RLVTT &               zx_dq_s_dh = dqsatl(tsurf(i), zx_qs)/RLVTT &
1599                    &               / zx_pkh(i)                    &               / zx_pkh(i)
1600            ENDIF            ENDIF
1601         ENDIF         ENDIF
# Line 1651  CONTAINS Line 1607  CONTAINS
1607      ENDDO      ENDDO
1608    
1609      ! === Calcul de la temperature de surface ===      ! === Calcul de la temperature de surface ===
1610      !  
1611      ! zx_sl = chaleur latente d'evaporation ou de sublimation      ! zx_sl = chaleur latente d'evaporation ou de sublimation
1612      !  
1613      do i = 1, knon      do i = 1, knon
1614         zx_sl(i) = RLVTT         zx_sl(i) = RLVTT
1615         if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT         if (tsurf(i) .LT. RTT) zx_sl(i) = RLSTT
# Line 1687  CONTAINS Line 1643  CONTAINS
1643      !  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
1644      bil_eau_s = (precip_rain * dtime) - (evap - snow_evap) * dtime      bil_eau_s = (precip_rain * dtime) - (evap - snow_evap) * dtime
1645    
1646      !  
1647      ! Y'a-t-il fonte de neige?      ! Y'a-t-il fonte de neige?
1648      !  
1649      ffonte=0.      ffonte=0.
1650      do i = 1, knon      do i = 1, knon
1651         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) &
1652              & .AND. tsurf_new(i) >= RTT)              & .AND. tsurf_new(i) >= RTT)
1653         if (neige_fond) then         if (neige_fond) then
1654            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))
1655            ffonte(i) = fq_fonte * RLMLT/dtime            ffonte(i) = fq_fonte * RLMLT/dtime
1656            snow(i) = max(0., snow(i) - fq_fonte)            snow(i) = max(0., snow(i) - fq_fonte)
1657            bil_eau_s(i) = bil_eau_s(i) + fq_fonte            bil_eau_s(i) = bil_eau_s(i) + fq_fonte
# Line 1703  CONTAINS Line 1659  CONTAINS
1659            !IM cf JLD OK                !IM cf JLD OK    
1660            !IM cf JLD/ GKtest fonte aussi pour la glace            !IM cf JLD/ GKtest fonte aussi pour la glace
1661            IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN            IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
1662               fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,0.0)               fq_fonte = MAX((tsurf_new(i)-RTT )/chaice, 0.0)
1663               ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime               ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime
1664               bil_eau_s(i) = bil_eau_s(i) + fq_fonte               bil_eau_s(i) = bil_eau_s(i) + fq_fonte
1665               tsurf_new(i) = RTT               tsurf_new(i) = RTT
1666            ENDIF            ENDIF
1667            d_ts(i) = tsurf_new(i) - tsurf(i)            d_ts(i) = tsurf_new(i) - tsurf(i)
1668         endif         endif
1669         !  
1670         !   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
1671         fqcalving(i) = max(0., snow(i) - snow_max)/dtime         fqcalving(i) = max(0., snow(i) - snow_max)/dtime
1672         snow(i)=min(snow(i),snow_max)         snow(i)=min(snow(i), snow_max)
1673         !  
1674         IF (nisurf == is_ter) then         IF (nisurf == is_ter) then
1675            qsol(i) = qsol(i) + bil_eau_s(i)            qsol(i) = qsol(i) + bil_eau_s(i)
1676            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)
Legend:
Removed from v.12  
changed lines
  Added in v.14
  ViewVC Help
Powered by ViewVC 1.1.21