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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 12 - (hide annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/interface_surf.f90
File size: 64594 byte(s) 
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

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