/[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 62 - (hide annotations)
Thu Jul 26 14:37:37 2012 UTC (11 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/Interface_surf/interfsurf_hq.f90
File size: 25989 byte(s) 
Changed handling of compiler in compilation system.

Removed the prefix letters "y", "p", "t" or "z" in some names of variables.

Replaced calls to NetCDF by calls to NetCDF95.

Extracted "ioget_calendar" procedures from "calendar.f90" into a
separate file.

Extracted to a separate file, "mathop2.f90", procedures that were not
part of the generic interface "mathop" in "mathop.f90".

Removed computation of "dq" in "bilan_dyn", which was not used.

In "iniadvtrac", removed schemes 20 Slopes and 30 Prather. Was not
compatible with declarations of array sizes.

In "clcdrag", "ustarhb", "vdif_kcay", "yamada4" and "coefkz", changed
the size of some arrays from "klon" to "knon".

Removed possible call to "conema3" in "physiq".

Removed unused argument "cd" in "yamada".
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 guez 62 real, intent(IN):: pctsrf(klon, nbsrf)
104 guez 54 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