1 |
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, intent(IN):: pctsrf(klon, nbsrf) |
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 |