/[lmdze]/trunk/Sources/phylmd/Interface_surf/interfsurf_hq.f
ViewVC logotype

Annotation of /trunk/Sources/phylmd/Interface_surf/interfsurf_hq.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 54 - (hide annotations)
Tue Dec 6 15:07:04 2011 UTC (12 years, 5 months ago) by guez
Original Path: trunk/libf/phylmd/Interface_surf/interfsurf_hq.f90
File size: 26001 byte(s) 
Removed Numerical Recipes procedure "ran1". Replaced calls to "ran1"
in "inidissip" by calls to intrinsic procedures.

Split file "interface_surf.f90" into a file with a module containing
only variables, "interface_surf", and single-procedure files. Gathered
files into directory "Interface_surf".

Added argument "cdivu" to "gradiv" and "gradiv2", "cdivh" to
"divgrad2" and "divgrad", and "crot" to "nxgraro2" and
"nxgrarot". "dissip" now uses variables "cdivu", "cdivh" and "crot"
from module "inidissip_m", so it can pass them to "gradiv2",
etc. Thanks to this modification, we avoid a circular dependency
betwwen "inidissip.f90" and "gradiv2.f90", etc. The value -1. used by
"gradiv2", for instance, during computation of eigenvalues is not the
value "cdivu" computed by "inidissip".

Extracted procedure "start_inter_3d" from module "startdyn", to its
own module.

In "inidissip", unrolled loop on "ii". I find it clearer now.

Moved variables "matriceun", "matriceus", "matricevn", "matricevs",
"matrinvn" and "matrinvs" from module "parafilt" to module
"inifilr_m". Moved variables "jfiltnu", "jfiltnv", "jfiltsu",
"jfiltsv" from module "coefils" to module "inifilr_m".
1 guez 54 module interfsurf_hq_m
2    
3     implicit none
4    
5     contains
6    
7     SUBROUTINE interfsurf_hq(itime, dtime, date0, jour, rmu0, klon, iim, jjm, &
8     nisurf, knon, knindex, pctsrf, rlon, rlat, cufi, cvfi, debut, lafin, &
9     ok_veget, soil_model, nsoilmx, tsoil, qsol, zlev, u1_lay, v1_lay, &
10     temp_air, spechum, epot_air, ccanopy, tq_cdrag, petAcoef, peqAcoef, &
11     petBcoef, peqBcoef, precip_rain, precip_snow, sollw, sollwdown, swnet, &
12     swdown, fder, taux, tauy, windsp, rugos, rugoro, albedo, snow, qsurf, &
13     tsurf, p1lay, ps, radsol, ocean, npas, nexca, zmasq, evap, fluxsens, &
14     fluxlat, dflux_l, dflux_s, tsol_rad, tsurf_new, alb_new, alblw, &
15     emis_new, z0_new, pctsrf_new, agesno, fqcalving, ffonte, &
16     run_off_lic_0, flux_o, flux_g, tslab, seaice)
17    
18     ! Cette routine sert d'aiguillage entre l'atmosphère et la surface
19     ! en général (sols continentaux, océans, glaces) pour les flux de
20     ! chaleur et d'humidité.
21     ! En pratique l'interface se fait entre la couche limite du modèle
22     ! atmosphérique ("clmain.F") et les routines de surface
23     ! ("sechiba", "oasis"...).
24    
25     ! L.Fairhead 02/2000
26    
27     use abort_gcm_m, only: abort_gcm
28     use gath_cpl, only: gath2cpl
29     use indicesol
30     use SUPHEC_M
31     use albsno_m, only: albsno
32     use interface_surf
33     use interfsur_lim_m, only: interfsur_lim
34     use calcul_fluxs_m, only: calcul_fluxs
35     use fonte_neige_m, only: fonte_neige
36     use interfoce_lim_m, only: interfoce_lim
37     use interfoce_slab_m, only: interfoce_slab
38    
39     ! Parametres d'entree
40     ! input:
41     ! klon nombre total de points de grille
42     ! iim, jjm nbres de pts de grille
43     ! dtime pas de temps de la physique (en s)
44     ! date0 jour initial
45     ! jour jour dans l'annee en cours,
46     ! rmu0 cosinus de l'angle solaire zenithal
47     ! nexca pas de temps couplage
48     ! nisurf index de la surface a traiter (1 = sol continental)
49     ! knon nombre de points de la surface a traiter
50     ! knindex index des points de la surface a traiter
51     ! pctsrf tableau des pourcentages de surface de chaque maille
52     ! rlon longitudes
53     ! rlat latitudes
54     ! cufi, cvfi resolution des mailles en x et y (m)
55     ! debut logical: 1er appel a la physique
56     ! lafin logical: dernier appel a la physique
57     ! ok_veget logical: appel ou non au schema de surface continental
58     ! (si false calcul simplifie des fluxs sur les continents)
59     ! zlev hauteur de la premiere couche
60     ! u1_lay vitesse u 1ere couche
61     ! v1_lay vitesse v 1ere couche
62     ! temp_air temperature de l'air 1ere couche
63     ! spechum humidite specifique 1ere couche
64     ! epot_air temp potentielle de l'air
65     ! ccanopy concentration CO2 canopee
66     ! tq_cdrag cdrag
67     ! petAcoef coeff. A de la resolution de la CL pour t
68     ! peqAcoef coeff. A de la resolution de la CL pour q
69     ! petBcoef coeff. B de la resolution de la CL pour t
70     ! peqBcoef coeff. B de la resolution de la CL pour q
71     ! precip_rain precipitation liquide
72     ! precip_snow precipitation solide
73     ! sollw flux IR net a la surface
74     ! sollwdown flux IR descendant a la surface
75     ! swnet flux solaire net
76     ! swdown flux solaire entrant a la surface
77     ! albedo albedo de la surface
78     ! tsurf temperature de surface
79     ! tslab temperature slab ocean
80     ! pctsrf_slab pourcentages (0-1) des sous-surfaces dans le slab
81     ! tmp_pctsrf_slab = pctsrf_slab
82     ! p1lay pression 1er niveau (milieu de couche)
83     ! ps pression au sol
84     ! radsol rayonnement net aus sol (LW + SW)
85     ! ocean type d'ocean utilise ("force" ou "slab" mais pas "couple")
86     ! fder derivee des flux (pour le couplage)
87     ! taux, tauy tension de vents
88     ! windsp module du vent a 10m
89     ! rugos rugosite
90     ! zmasq masque terre/ocean
91     ! rugoro rugosite orographique
92     ! run_off_lic_0 runoff glacier du pas de temps precedent
93     integer, intent(IN) :: itime ! numero du pas de temps
94     integer, intent(IN) :: iim, jjm
95     integer, intent(IN) :: klon
96     real, intent(IN) :: dtime
97     real, intent(IN) :: date0
98     integer, intent(IN) :: jour
99     real, intent(IN) :: rmu0(klon)
100     integer, intent(IN) :: nisurf
101     integer, intent(IN) :: knon
102     integer, dimension(klon), intent(in) :: knindex
103     real, dimension(klon, nbsrf), intent(IN) :: pctsrf
104     logical, intent(IN) :: debut, lafin, ok_veget
105     real, dimension(klon), intent(IN) :: rlon, rlat
106     real, dimension(klon), intent(IN) :: cufi, cvfi
107     real, dimension(klon), intent(INOUT) :: tq_cdrag
108     real, dimension(klon), intent(IN) :: zlev
109     real, dimension(klon), intent(IN) :: u1_lay, v1_lay
110     real, dimension(klon), intent(IN) :: temp_air, spechum
111     real, dimension(klon), intent(IN) :: epot_air, ccanopy
112     real, dimension(klon), intent(IN) :: petAcoef, peqAcoef
113     real, dimension(klon), intent(IN) :: petBcoef, peqBcoef
114     real, dimension(klon), intent(IN) :: precip_rain, precip_snow
115     real, dimension(klon), intent(IN) :: sollw, sollwdown, swnet, swdown
116     real, dimension(klon), intent(IN) :: ps, albedo
117     real, dimension(klon), intent(IN) :: tsurf, p1lay
118     !IM: "slab" ocean
119     real, dimension(klon), intent(INOUT) :: tslab
120     real, allocatable, dimension(:), save :: tmp_tslab
121     real, dimension(klon), intent(OUT) :: flux_o, flux_g
122     real, dimension(klon), intent(INOUT) :: seaice ! glace de mer (kg/m2)
123     REAL, DIMENSION(klon), INTENT(INOUT) :: radsol, fder
124     real, dimension(klon), intent(IN) :: zmasq
125     real, dimension(klon), intent(IN) :: taux, tauy, rugos, rugoro
126     real, dimension(klon), intent(IN) :: windsp
127     character(len=*), intent(IN):: ocean
128     integer :: npas, nexca ! nombre et pas de temps couplage
129     real, dimension(klon), intent(INOUT) :: evap, snow, qsurf
130     !! PB ajout pour soil
131     logical, intent(in):: soil_model
132     integer :: nsoilmx
133     REAL, DIMENSION(klon, nsoilmx) :: tsoil
134     REAL, dimension(klon), intent(INOUT) :: qsol
135     REAL, dimension(klon) :: soilcap
136     REAL, dimension(klon) :: soilflux
137    
138     ! Parametres de sortie
139     ! output:
140     ! evap evaporation totale
141     ! fluxsens flux de chaleur sensible
142     ! fluxlat flux de chaleur latente
143     ! tsol_rad
144     ! tsurf_new temperature au sol
145     ! alb_new albedo
146     ! emis_new emissivite
147     ! z0_new surface roughness
148     ! pctsrf_new nouvelle repartition des surfaces
149     real, dimension(klon), intent(OUT):: fluxsens, fluxlat
150     real, dimension(klon), intent(OUT):: tsol_rad, tsurf_new, alb_new
151     real, dimension(klon), intent(OUT):: alblw
152     real, dimension(klon), intent(OUT):: emis_new, z0_new
153     real, dimension(klon), intent(OUT):: dflux_l, dflux_s
154     real, dimension(klon, nbsrf), intent(OUT) :: pctsrf_new
155     real, dimension(klon), intent(INOUT):: agesno
156     real, dimension(klon), intent(INOUT):: run_off_lic_0
157    
158     ! Flux thermique utiliser pour fondre la neige
159     !jld a rajouter real, dimension(klon), intent(INOUT):: ffonte
160     real, dimension(klon), intent(INOUT):: ffonte
161     ! Flux d'eau "perdue" par la surface et nécessaire pour que limiter la
162     ! hauteur de neige, en kg/m2/s
163     !jld a rajouter real, dimension(klon), intent(INOUT):: fqcalving
164     real, dimension(klon), intent(INOUT):: fqcalving
165     !IM: "slab" ocean - Local
166     real, parameter :: t_grnd=271.35
167     real, dimension(klon) :: zx_sl
168     integer i
169     real, allocatable, dimension(:), save :: tmp_flux_o, tmp_flux_g
170     real, allocatable, dimension(:), save :: tmp_radsol
171     real, allocatable, dimension(:, :), save :: tmp_pctsrf_slab
172     real, allocatable, dimension(:), save :: tmp_seaice
173    
174     ! Local
175     character (len = 20), save :: modname = 'interfsurf_hq'
176     character (len = 80) :: abort_message
177     logical, save :: first_call = .true.
178     integer, save :: error
179     integer :: ii
180     logical, save :: check = .false.
181     real, dimension(klon):: cal, beta, dif_grnd, capsol
182     real, parameter :: calice=1.0/(5.1444e+06*0.15), tau_gl=86400.*5.
183     real, parameter :: calsno=1./(2.3867e+06*.15)
184     real, dimension(klon):: tsurf_temp
185     real, dimension(klon):: alb_neig, alb_eau
186     real, DIMENSION(klon):: zfra
187     logical :: cumul = .false.
188     INTEGER, dimension(1) :: iloc
189     real, dimension(klon):: fder_prev
190     REAL, dimension(klon) :: bidule
191    
192     !-------------------------------------------------------------
193    
194     if (check) write(*, *) 'Entree ', modname
195    
196     ! On doit commencer par appeler les schemas de surfaces continentales
197     ! car l'ocean a besoin du ruissellement qui est y calcule
198    
199     if (first_call) then
200     call conf_interface(tau_calv)
201     if (nisurf /= is_ter .and. klon > 1) then
202     write(*, *)' *** Warning ***'
203     write(*, *)' nisurf = ', nisurf, ' /= is_ter = ', is_ter
204     write(*, *)'or on doit commencer par les surfaces continentales'
205     abort_message='voir ci-dessus'
206     call abort_gcm(modname, abort_message, 1)
207     endif
208     if (ocean /= 'slab' .and. ocean /= 'force') then
209     write(*, *)' *** Warning ***'
210     write(*, *)'Option couplage pour l''ocean = ', ocean
211     abort_message='option pour l''ocean non valable'
212     call abort_gcm(modname, abort_message, 1)
213     endif
214     if ( is_oce > is_sic ) then
215     write(*, *)' *** Warning ***'
216     write(*, *)' Pour des raisons de sequencement dans le code'
217     write(*, *)' l''ocean doit etre traite avant la banquise'
218     write(*, *)' or is_oce = ', is_oce, '> is_sic = ', is_sic
219     abort_message='voir ci-dessus'
220     call abort_gcm(modname, abort_message, 1)
221     endif
222     endif
223     first_call = .false.
224    
225     ! Initialisations diverses
226    
227     ffonte(1:knon)=0.
228     fqcalving(1:knon)=0.
229    
230     cal = 999999.
231     beta = 999999.
232     dif_grnd = 999999.
233     capsol = 999999.
234     alb_new = 999999.
235     z0_new = 999999.
236     alb_neig = 999999.
237     tsurf_new = 999999.
238     alblw = 999999.
239    
240     !IM: "slab" ocean; initialisations
241     flux_o = 0.
242     flux_g = 0.
243    
244     if (.not. allocated(tmp_flux_o)) then
245     allocate(tmp_flux_o(klon), stat = error)
246     DO i=1, knon
247     tmp_flux_o(knindex(i))=flux_o(i)
248     ENDDO
249     if (error /= 0) then
250     abort_message='Pb allocation tmp_flux_o'
251     call abort_gcm(modname, abort_message, 1)
252     endif
253     endif
254     if (.not. allocated(tmp_flux_g)) then
255     allocate(tmp_flux_g(klon), stat = error)
256     DO i=1, knon
257     tmp_flux_g(knindex(i))=flux_g(i)
258     ENDDO
259     if (error /= 0) then
260     abort_message='Pb allocation tmp_flux_g'
261     call abort_gcm(modname, abort_message, 1)
262     endif
263     endif
264     if (.not. allocated(tmp_radsol)) then
265     allocate(tmp_radsol(klon), stat = error)
266     if (error /= 0) then
267     abort_message='Pb allocation tmp_radsol'
268     call abort_gcm(modname, abort_message, 1)
269     endif
270     endif
271     DO i=1, knon
272     tmp_radsol(knindex(i))=radsol(i)
273     ENDDO
274     if (.not. allocated(tmp_pctsrf_slab)) then
275     allocate(tmp_pctsrf_slab(klon, nbsrf), stat = error)
276     if (error /= 0) then
277     abort_message='Pb allocation tmp_pctsrf_slab'
278     call abort_gcm(modname, abort_message, 1)
279     endif
280     DO i=1, klon
281     tmp_pctsrf_slab(i, 1:nbsrf)=pctsrf(i, 1:nbsrf)
282     ENDDO
283     endif
284    
285     if (.not. allocated(tmp_seaice)) then
286     allocate(tmp_seaice(klon), stat = error)
287     if (error /= 0) then
288     abort_message='Pb allocation tmp_seaice'
289     call abort_gcm(modname, abort_message, 1)
290     endif
291     DO i=1, klon
292     tmp_seaice(i)=seaice(i)
293     ENDDO
294     endif
295    
296     if (.not. allocated(tmp_tslab)) then
297     allocate(tmp_tslab(klon), stat = error)
298     if (error /= 0) then
299     abort_message='Pb allocation tmp_tslab'
300     call abort_gcm(modname, abort_message, 1)
301     endif
302     endif
303     DO i=1, klon
304     tmp_tslab(i)=tslab(i)
305     ENDDO
306    
307     ! Aiguillage vers les differents schemas de surface
308    
309     if (nisurf == is_ter) then
310    
311     ! Surface "terre" appel a l'interface avec les sols continentaux
312    
313     ! allocation du run-off
314     if (.not. allocated(coastalflow)) then
315     allocate(coastalflow(knon), stat = error)
316     if (error /= 0) then
317     abort_message='Pb allocation coastalflow'
318     call abort_gcm(modname, abort_message, 1)
319     endif
320     allocate(riverflow(knon), stat = error)
321     if (error /= 0) then
322     abort_message='Pb allocation riverflow'
323     call abort_gcm(modname, abort_message, 1)
324     endif
325     allocate(run_off(knon), stat = error)
326     if (error /= 0) then
327     abort_message='Pb allocation run_off'
328     call abort_gcm(modname, abort_message, 1)
329     endif
330     !cym
331     run_off=0.0
332     !cym
333    
334     !!$PB
335     ALLOCATE (tmp_rriv(iim, jjm+1), stat=error)
336     if (error /= 0) then
337     abort_message='Pb allocation tmp_rriv'
338     call abort_gcm(modname, abort_message, 1)
339     endif
340     ALLOCATE (tmp_rcoa(iim, jjm+1), stat=error)
341     if (error /= 0) then
342     abort_message='Pb allocation tmp_rcoa'
343     call abort_gcm(modname, abort_message, 1)
344     endif
345     ALLOCATE (tmp_rlic(iim, jjm+1), stat=error)
346     if (error /= 0) then
347     abort_message='Pb allocation tmp_rlic'
348     call abort_gcm(modname, abort_message, 1)
349     endif
350     tmp_rriv = 0.0
351     tmp_rcoa = 0.0
352     tmp_rlic = 0.0
353    
354     !!$
355     else if (size(coastalflow) /= knon) then
356     write(*, *)'Bizarre, le nombre de points continentaux'
357     write(*, *)'a change entre deux appels. J''arrete ...'
358     abort_message='voir ci-dessus'
359     call abort_gcm(modname, abort_message, 1)
360     endif
361     coastalflow = 0.
362     riverflow = 0.
363    
364     ! Calcul age de la neige
365    
366     if (.not. ok_veget) then
367     ! calcul albedo: lecture albedo fichier boundary conditions
368     ! puis ajout albedo neige
369     call interfsur_lim(itime, dtime, jour, klon, nisurf, knon, knindex, &
370     debut, alb_new, z0_new)
371    
372     ! calcul snow et qsurf, hydrol adapté
373     CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
374    
375     IF (soil_model) THEN
376     CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, &
377     soilflux)
378     cal(1:knon) = RCPD / soilcap(1:knon)
379     radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
380     ELSE
381     cal = RCPD * capsol
382     ENDIF
383     CALL calcul_fluxs( klon, knon, nisurf, dtime, &
384     tsurf, p1lay, cal, beta, tq_cdrag, ps, &
385     precip_rain, precip_snow, snow, qsurf, &
386     radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
387     petAcoef, peqAcoef, petBcoef, peqBcoef, &
388     tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
389    
390     CALL fonte_neige( klon, knon, nisurf, dtime, &
391     tsurf, p1lay, cal, beta, tq_cdrag, ps, &
392     precip_rain, precip_snow, snow, qsol, &
393     radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
394     petAcoef, peqAcoef, petBcoef, peqBcoef, &
395     tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
396     fqcalving, ffonte, run_off_lic_0)
397    
398     call albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)
399     where (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
400     zfra(1:knon) = max(0.0, min(1.0, snow(1:knon)/(snow(1:knon)+10.0)))
401     alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
402     alb_new(1 : knon)*(1.0-zfra(1:knon))
403     z0_new = sqrt(z0_new**2+rugoro**2)
404     alblw(1 : knon) = alb_new(1 : knon)
405     endif
406    
407     ! Remplissage des pourcentages de surface
408     pctsrf_new(:, nisurf) = pctsrf(:, nisurf)
409     else if (nisurf == is_oce) then
410     ! Surface "ocean" appel a l'interface avec l'ocean
411     if (ocean == 'slab') then
412     tsurf_new = tsurf
413     pctsrf_new = tmp_pctsrf_slab
414     else
415     ! lecture conditions limites
416     call interfoce_lim(itime, dtime, jour, klon, nisurf, knon, knindex, &
417     debut, tsurf_new, pctsrf_new)
418     endif
419    
420     tsurf_temp = tsurf_new
421     cal = 0.
422     beta = 1.
423     dif_grnd = 0.
424     alb_neig = 0.
425     agesno = 0.
426    
427     call calcul_fluxs( klon, knon, nisurf, dtime, &
428     tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
429     precip_rain, precip_snow, snow, qsurf, &
430     radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
431     petAcoef, peqAcoef, petBcoef, peqBcoef, &
432     tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
433    
434     fder_prev = fder
435     fder = fder_prev + dflux_s + dflux_l
436    
437     iloc = maxloc(fder(1:klon))
438     if (check.and.fder(iloc(1))> 0.) then
439     WRITE(*, *)'**** Debug fder****'
440     WRITE(*, *)'max fder(', iloc(1), ') = ', fder(iloc(1))
441     WRITE(*, *)'fder_prev, dflux_s, dflux_l', fder_prev(iloc(1)), &
442     dflux_s(iloc(1)), dflux_l(iloc(1))
443     endif
444    
445     !IM: flux ocean-atmosphere utile pour le "slab" ocean
446     DO i=1, knon
447     zx_sl(i) = RLVTT
448     if (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
449     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
450     tmp_flux_o(knindex(i)) = flux_o(i)
451     tmp_radsol(knindex(i))=radsol(i)
452     ENDDO
453    
454     ! 2eme appel a interfoce pour le cumul des champs (en particulier
455     ! fluxsens et fluxlat calcules dans calcul_fluxs)
456    
457     if (ocean == 'slab ') then
458     seaice=tmp_seaice
459     cumul = .true.
460     call interfoce_slab(klon, debut, itime, dtime, jour, &
461     tmp_radsol, tmp_flux_o, tmp_flux_g, pctsrf, &
462     tslab, seaice, pctsrf_new)
463    
464     tmp_pctsrf_slab=pctsrf_new
465     DO i=1, knon
466     tsurf_new(i)=tslab(knindex(i))
467     ENDDO
468     endif
469    
470     ! calcul albedo
471     if ( minval(rmu0) == maxval(rmu0) .and. minval(rmu0) == -999.999 ) then
472     CALL alboc(FLOAT(jour), rlat, alb_eau)
473     else ! cycle diurne
474     CALL alboc_cd(rmu0, alb_eau)
475     endif
476     DO ii =1, knon
477     alb_new(ii) = alb_eau(knindex(ii))
478     enddo
479    
480     z0_new = sqrt(rugos**2 + rugoro**2)
481     alblw(1:knon) = alb_new(1:knon)
482     else if (nisurf == is_sic) then
483     if (check) write(*, *)'sea ice, nisurf = ', nisurf
484    
485     ! Surface "glace de mer" appel a l'interface avec l'ocean
486    
487    
488     if (ocean == 'slab ') then
489     pctsrf_new=tmp_pctsrf_slab
490    
491     DO ii = 1, knon
492     tsurf_new(ii) = tsurf(ii)
493     IF (pctsrf_new(knindex(ii), nisurf) < EPSFRA) then
494     snow(ii) = 0.0
495     tsurf_new(ii) = RTT - 1.8
496     IF (soil_model) tsoil(ii, :) = RTT -1.8
497     ENDIF
498     ENDDO
499    
500     CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
501    
502     IF (soil_model) THEN
503     CALL soil(dtime, nisurf, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
504     cal(1:knon) = RCPD / soilcap(1:knon)
505     radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
506     ELSE
507     dif_grnd = 1.0 / tau_gl
508     cal = RCPD * calice
509     WHERE (snow > 0.0) cal = RCPD * calsno
510     ENDIF
511     tsurf_temp = tsurf_new
512     beta = 1.0
513    
514     ELSE
515     ! ! lecture conditions limites
516     CALL interfoce_lim(itime, dtime, jour, &
517     klon, nisurf, knon, knindex, &
518     debut, &
519     tsurf_new, pctsrf_new)
520    
521     !IM cf LF
522     DO ii = 1, knon
523     tsurf_new(ii) = tsurf(ii)
524     !IMbad IF (pctsrf_new(ii, nisurf) < EPSFRA) then
525     IF (pctsrf_new(knindex(ii), nisurf) < EPSFRA) then
526     snow(ii) = 0.0
527     !IM cf LF/JLD tsurf(ii) = RTT - 1.8
528     tsurf_new(ii) = RTT - 1.8
529     IF (soil_model) tsoil(ii, :) = RTT -1.8
530     endif
531     enddo
532    
533     CALL calbeta(dtime, nisurf, knon, snow, qsol, beta, capsol, dif_grnd)
534    
535     IF (soil_model) THEN
536     !IM cf LF/JLD CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux)
537     CALL soil(dtime, nisurf, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
538     cal(1:knon) = RCPD / soilcap(1:knon)
539     radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
540     dif_grnd = 0.
541     ELSE
542     dif_grnd = 1.0 / tau_gl
543     cal = RCPD * calice
544     WHERE (snow > 0.0) cal = RCPD * calsno
545     ENDIF
546     !IMbadtsurf_temp = tsurf
547     tsurf_temp = tsurf_new
548     beta = 1.0
549     ENDIF
550    
551     CALL calcul_fluxs( klon, knon, nisurf, dtime, &
552     tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
553     precip_rain, precip_snow, snow, qsurf, &
554     radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
555     petAcoef, peqAcoef, petBcoef, peqBcoef, &
556     tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
557    
558     !IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
559     DO i = 1, knon
560     flux_g(i) = 0.0
561    
562     !IM: faire dependre le coefficient de conduction de la glace de mer
563     ! de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
564     ! actuel correspond a 3m de glace de mer, cf. L.Li
565    
566     ! IF(1.EQ.0) THEN
567     ! IF(siceh(i).GT.0.) THEN
568     ! new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)
569     ! ELSE
570     ! new_dif_grnd(i) = 0.
571     ! ENDIF
572     ! ENDIF !(1.EQ.0) THEN
573    
574     IF (cal(i).GT.1.0e-15) flux_g(i)=(tsurf_new(i)-t_grnd) &
575     * dif_grnd(i) *RCPD/cal(i)
576     ! & * new_dif_grnd(i) *RCPD/cal(i)
577     tmp_flux_g(knindex(i))=flux_g(i)
578     tmp_radsol(knindex(i))=radsol(i)
579     ENDDO
580    
581     CALL fonte_neige( klon, knon, nisurf, dtime, &
582     tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
583     precip_rain, precip_snow, snow, qsol, &
584     radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
585     petAcoef, peqAcoef, petBcoef, peqBcoef, &
586     tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
587     fqcalving, ffonte, run_off_lic_0)
588    
589     ! calcul albedo
590    
591     CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)
592     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
593     zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon)/(snow(1:knon)+10.0)))
594     alb_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
595     0.6 * (1.0-zfra(1:knon))
596    
597     fder_prev = fder
598     fder = fder_prev + dflux_s + dflux_l
599    
600     iloc = maxloc(fder(1:klon))
601     if (check.and.fder(iloc(1))> 0.) then
602     WRITE(*, *)'**** Debug fder ****'
603     WRITE(*, *)'max fder(', iloc(1), ') = ', fder(iloc(1))
604     WRITE(*, *)'fder_prev, dflux_s, dflux_l', fder_prev(iloc(1)), &
605     dflux_s(iloc(1)), dflux_l(iloc(1))
606     endif
607    
608    
609     ! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
610    
611     z0_new = 0.002
612     z0_new = SQRT(z0_new**2+rugoro**2)
613     alblw(1:knon) = alb_new(1:knon)
614    
615     else if (nisurf == is_lic) then
616    
617     if (check) write(*, *)'glacier, nisurf = ', nisurf
618    
619     if (.not. allocated(run_off_lic)) then
620     allocate(run_off_lic(knon), stat = error)
621     if (error /= 0) then
622     abort_message='Pb allocation run_off_lic'
623     call abort_gcm(modname, abort_message, 1)
624     endif
625     run_off_lic = 0.
626     endif
627    
628     ! Surface "glacier continentaux" appel a l'interface avec le sol
629    
630     IF (soil_model) THEN
631     CALL soil(dtime, nisurf, knon, snow, tsurf, tsoil, soilcap, soilflux)
632     cal(1:knon) = RCPD / soilcap(1:knon)
633     radsol(1:knon) = radsol(1:knon) + soilflux(1:knon)
634     ELSE
635     cal = RCPD * calice
636     WHERE (snow > 0.0) cal = RCPD * calsno
637     ENDIF
638     beta = 1.0
639     dif_grnd = 0.0
640    
641     call calcul_fluxs( klon, knon, nisurf, dtime, &
642     tsurf, p1lay, cal, beta, tq_cdrag, ps, &
643     precip_rain, precip_snow, snow, qsurf, &
644     radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
645     petAcoef, peqAcoef, petBcoef, peqBcoef, &
646     tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
647    
648     call fonte_neige( klon, knon, nisurf, dtime, &
649     tsurf, p1lay, cal, beta, tq_cdrag, ps, &
650     precip_rain, precip_snow, snow, qsol, &
651     radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
652     petAcoef, peqAcoef, petBcoef, peqBcoef, &
653     tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
654     fqcalving, ffonte, run_off_lic_0)
655    
656     ! passage du run-off des glaciers calcule dans fonte_neige au coupleur
657     bidule=0.
658     bidule(1:knon)= run_off_lic(1:knon)
659     call gath2cpl(bidule, tmp_rlic, klon, knon, iim, jjm, knindex)
660    
661     ! calcul albedo
662    
663     CALL albsno(klon, knon, dtime, agesno, alb_neig, precip_snow)
664     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
665     zfra(1:knon) = MAX(0.0, MIN(1.0, snow(1:knon)/(snow(1:knon)+10.0)))
666     alb_new(1 : knon) = alb_neig(1 : knon)*zfra(1:knon) + &
667     0.6 * (1.0-zfra(1:knon))
668    
669     !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux"
670     ! alb_new(1 : knon) = 0.6 !IM cf FH/GK
671     ! alb_new(1 : knon) = 0.82
672     ! alb_new(1 : knon) = 0.77 !211003 Ksta0.77
673     ! alb_new(1 : knon) = 0.8 !KstaTER0.8 & LMD_ARMIP5
674     !IM: KstaTER0.77 & LMD_ARMIP6
675     alb_new(1 : knon) = 0.77
676    
677    
678     ! Rugosite
679    
680     z0_new = rugoro
681    
682     ! Remplissage des pourcentages de surface
683    
684     pctsrf_new(:, nisurf) = pctsrf(:, nisurf)
685    
686     alblw(1:knon) = alb_new(1:knon)
687     else
688     write(*, *)'Index surface = ', nisurf
689     abort_message = 'Index surface non valable'
690     call abort_gcm(modname, abort_message, 1)
691     endif
692    
693     END SUBROUTINE interfsurf_hq
694    
695     end module interfsurf_hq_m
  ViewVC Help
Powered by ViewVC 1.1.21