source: LMDZ5/trunk/libf/phylmd/fonte_neige_mod.F90 @ 1635

Last change on this file since 1635 was 1504, checked in by Laurent Fairhead, 14 years ago

Bug in runoff caclulation over land ice.

Bug sur le calcul du ruisellement sur les glaciers
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 12.3 KB
Line 
1!
2! $Header$
3!
4MODULE fonte_neige_mod
5!
6! This module will treat the process of snow, melting, accumulating, calving, in
7! case of simplified soil model.
8!
9!****************************************************************************************
10  USE dimphy, ONLY : klon
11
12  IMPLICIT NONE
13  SAVE
14
15! run_off_ter and run_off_lic are the runoff at the compressed grid knon for
16! land and land-ice respectively
17! Note: run_off_lic is used in mod_landice and therfore not private
18  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE    :: run_off_ter
19  !$OMP THREADPRIVATE(run_off_ter)
20  REAL, ALLOCATABLE, DIMENSION(:)             :: run_off_lic
21  !$OMP THREADPRIVATE(run_off_lic)
22
23! run_off_lic_0 is the runoff at land-ice a time-step earlier, on the global 1D array grid
24  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE    :: run_off_lic_0
25  !$OMP THREADPRIVATE(run_off_lic_0)
26 
27  REAL, PRIVATE                               :: tau_calv 
28  !$OMP THREADPRIVATE(tau_calv)
29  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: ffonte_global
30  !$OMP THREADPRIVATE(ffonte_global)
31  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: fqfonte_global
32  !$OMP THREADPRIVATE(fqfonte_global)
33  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE  :: fqcalving_global
34  !$OMP THREADPRIVATE(fqcalving_global)
35
36CONTAINS
37!
38!****************************************************************************************
39!
40  SUBROUTINE fonte_neige_init(restart_runoff)
41
42! This subroutine allocates and initialize variables in the module.
43! The variable run_off_lic_0 is initialized to the field read from
44! restart file. The other variables are initialized to zero.
45!
46    INCLUDE "indicesol.h"
47!****************************************************************************************
48! Input argument
49    REAL, DIMENSION(klon), INTENT(IN) :: restart_runoff
50
51! Local variables
52    INTEGER                           :: error
53    CHARACTER (len = 80)              :: abort_message
54    CHARACTER (len = 20)              :: modname = 'fonte_neige_init'
55
56
57!****************************************************************************************
58! Allocate run-off at landice and initilize with field read from restart
59!
60!****************************************************************************************
61
62    ALLOCATE(run_off_lic_0(klon), stat = error)
63    IF (error /= 0) THEN
64       abort_message='Pb allocation run_off_lic'
65       CALL abort_gcm(modname,abort_message,1)
66    ENDIF
67    run_off_lic_0(:) = restart_runoff(:)
68
69!****************************************************************************************
70! Allocate other variables and initilize to zero
71!
72!****************************************************************************************
73    ALLOCATE(run_off_ter(klon), stat = error)
74    IF (error /= 0) THEN
75       abort_message='Pb allocation run_off_ter'
76       CALL abort_gcm(modname,abort_message,1)
77    ENDIF
78    run_off_ter(:) = 0.
79   
80    ALLOCATE(run_off_lic(klon), stat = error)
81    IF (error /= 0) THEN
82       abort_message='Pb allocation run_off_lic'
83       CALL abort_gcm(modname,abort_message,1)
84    ENDIF
85    run_off_lic(:) = 0.
86   
87    ALLOCATE(ffonte_global(klon,nbsrf))
88    IF (error /= 0) THEN
89       abort_message='Pb allocation ffonte_global'
90       CALL abort_gcm(modname,abort_message,1)
91    ENDIF
92    ffonte_global(:,:) = 0.0
93
94    ALLOCATE(fqfonte_global(klon,nbsrf))
95    IF (error /= 0) THEN
96       abort_message='Pb allocation fqfonte_global'
97       CALL abort_gcm(modname,abort_message,1)
98    ENDIF
99    fqfonte_global(:,:) = 0.0
100
101    ALLOCATE(fqcalving_global(klon,nbsrf))
102    IF (error /= 0) THEN
103       abort_message='Pb allocation fqcalving_global'
104       CALL abort_gcm(modname,abort_message,1)
105    ENDIF
106    fqcalving_global(:,:) = 0.0
107
108!****************************************************************************************
109! Read tau_calv
110!
111!****************************************************************************************
112    CALL conf_interface(tau_calv)
113
114
115  END SUBROUTINE fonte_neige_init
116!
117!****************************************************************************************
118!
119  SUBROUTINE fonte_neige( knon, nisurf, knindex, dtime, &
120       tsurf, precip_rain, precip_snow, &
121       snow, qsol, tsurf_new, evap)
122       
123! Routine de traitement de la fonte de la neige dans le cas du traitement
124! de sol simplifie!
125! LF 03/2001
126! input:
127!   knon         nombre de points a traiter
128!   nisurf       surface a traiter
129!   knindex      index des mailles valables pour surface a traiter
130!   dtime       
131!   tsurf        temperature de surface
132!   precip_rain  precipitations liquides
133!   precip_snow  precipitations solides
134!
135! input/output:
136!   snow         champs hauteur de neige
137!   qsol         hauteur d'eau contenu dans le sol
138!   tsurf_new    temperature au sol
139!   evap
140!
141  INCLUDE "indicesol.h"
142  INCLUDE "dimensions.h"
143  INCLUDE "YOETHF.h"
144  INCLUDE "YOMCST.h"
145  INCLUDE "FCTTRE.h"
146  INCLUDE "clesphys.h"
147
148! Input variables
149!****************************************************************************************
150    INTEGER, INTENT(IN)                  :: knon
151    INTEGER, INTENT(IN)                  :: nisurf
152    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
153    REAL   , INTENT(IN)                  :: dtime
154    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf
155    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain
156    REAL, DIMENSION(klon), INTENT(IN)    :: precip_snow
157   
158! Input/Output variables
159!****************************************************************************************
160
161    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
162    REAL, DIMENSION(klon), INTENT(INOUT) :: qsol
163    REAL, DIMENSION(klon), INTENT(INOUT) :: tsurf_new
164    REAL, DIMENSION(klon), INTENT(INOUT) :: evap
165
166! Local variables
167!****************************************************************************************
168
169    INTEGER               :: i, j
170    REAL                  :: fq_fonte
171    REAL                  :: coeff_rel
172    REAL, PARAMETER       :: snow_max=3000.
173    REAL, PARAMETER       :: max_eau_sol = 150.0
174!! PB temporaire en attendant mieux pour le modele de neige
175! REAL, parameter :: chasno = RLMLT/(2.3867E+06*0.15)
176    REAL, PARAMETER       :: chasno = 3.334E+05/(2.3867E+06*0.15)
177!IM cf JLD/ GKtest
178    REAL, PARAMETER       :: chaice = 3.334E+05/(2.3867E+06*0.15)
179! fin GKtest
180    REAL, DIMENSION(klon) :: ffonte
181    REAL, DIMENSION(klon) :: fqcalving, fqfonte
182    REAL, DIMENSION(klon) :: d_ts
183    REAL, DIMENSION(klon) :: bil_eau_s, snow_evap
184
185    LOGICAL               :: neige_fond
186
187!****************************************************************************************
188! Start calculation
189! - Initialization
190!
191!****************************************************************************************
192    coeff_rel = dtime/(tau_calv * rday)
193   
194    bil_eau_s(:) = 0.
195
196!****************************************************************************************
197! - Increment snow due to precipitation and evaporation
198! - Calculate the water balance due to precipitation and evaporation (bil_eau_s)
199!
200!****************************************************************************************
201    WHERE (precip_snow > 0.)
202       snow = snow + (precip_snow * dtime)
203    END WHERE
204
205    snow_evap = 0.
206    WHERE (evap > 0. )
207       snow_evap = MIN (snow / dtime, evap)
208       snow = snow - snow_evap * dtime
209       snow = MAX(0.0, snow)
210    END WHERE
211   
212    bil_eau_s(:) = (precip_rain(:) * dtime) - (evap(:) - snow_evap(:)) * dtime
213
214
215!****************************************************************************************
216! - Calculate melting snow
217! - Calculate calving and decrement snow, if there are to much snow
218! - Update temperature at surface
219!
220!****************************************************************************************
221
222    ffonte(:) = 0.0
223    fqcalving(:) = 0.0
224    fqfonte(:) = 0.0
225    DO i = 1, knon
226       ! Y'a-t-il fonte de neige?
227       neige_fond = ((snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) &
228            .AND. tsurf_new(i) >= RTT)
229       IF (neige_fond) THEN
230          fq_fonte     = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i))
231          ffonte(i)    = fq_fonte * RLMLT/dtime
232          fqfonte(i)   = fq_fonte/dtime
233          snow(i)      = MAX(0., snow(i) - fq_fonte)
234          bil_eau_s(i) = bil_eau_s(i) + fq_fonte
235          tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno 
236
237!IM cf JLD OK     
238!IM cf JLD/ GKtest fonte aussi pour la glace
239          IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
240             fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,0.0)
241             ffonte(i) = ffonte(i) + fq_fonte * RLMLT/dtime
242             IF ( ok_lic_melt ) THEN
243                fqfonte(i) = fqfonte(i) + fq_fonte/dtime
244                bil_eau_s(i) = bil_eau_s(i) + fq_fonte
245             ENDIF
246             tsurf_new(i) = RTT
247          ENDIF
248          d_ts(i) = tsurf_new(i) - tsurf(i)
249       ENDIF
250
251       ! s'il y a une hauteur trop importante de neige, elle s'coule
252       fqcalving(i) = MAX(0., snow(i) - snow_max)/dtime
253       snow(i)=MIN(snow(i),snow_max)
254    END DO
255
256
257    IF (nisurf == is_ter) THEN
258       DO i = 1, knon
259          qsol(i) = qsol(i) + bil_eau_s(i)
260          run_off_ter(i) = run_off_ter(i) + MAX(qsol(i) - max_eau_sol, 0.0)
261          qsol(i) = MIN(qsol(i), max_eau_sol)
262       END DO
263    ELSE IF (nisurf == is_lic) THEN
264       DO i = 1, knon
265          j = knindex(i)
266          run_off_lic(i)   = (coeff_rel *  fqcalving(i)) + &
267               (1. - coeff_rel) * run_off_lic_0(j)
268          run_off_lic_0(j) = run_off_lic(i)
269          run_off_lic(i)   = run_off_lic(i) + fqfonte(i) + precip_rain(i)
270       END DO
271    ENDIF
272   
273!****************************************************************************************
274! Save ffonte, fqfonte and fqcalving in global arrays for each
275! sub-surface separately
276!
277!****************************************************************************************
278    DO i = 1, knon
279       ffonte_global(knindex(i),nisurf)    = ffonte(i)
280       fqfonte_global(knindex(i),nisurf)   = fqfonte(i)
281       fqcalving_global(knindex(i),nisurf) = fqcalving(i)
282    ENDDO
283
284  END SUBROUTINE fonte_neige
285!
286!****************************************************************************************
287!
288  SUBROUTINE fonte_neige_final(restart_runoff)
289!
290! This subroutine returns run_off_lic_0 for later writing to restart file.
291!
292!****************************************************************************************
293    REAL, DIMENSION(klon), INTENT(OUT) :: restart_runoff
294
295!****************************************************************************************
296! Set the output variables
297    restart_runoff(:) = run_off_lic_0(:)
298
299! Deallocation of all varaibles in the module
300!   DEALLOCATE(run_off_lic_0, run_off_ter, run_off_lic, ffonte_global, &
301!        fqfonte_global, fqcalving_global)
302
303    IF (ALLOCATED(run_off_lic_0)) DEALLOCATE(run_off_lic_0)
304    IF (ALLOCATED(run_off_ter)) DEALLOCATE(run_off_ter)
305    IF (ALLOCATED(run_off_lic)) DEALLOCATE(run_off_lic)
306    IF (ALLOCATED(ffonte_global)) DEALLOCATE(ffonte_global)
307    IF (ALLOCATED(fqfonte_global)) DEALLOCATE(fqfonte_global)
308    IF (ALLOCATED(fqcalving_global)) DEALLOCATE(fqcalving_global)
309
310  END SUBROUTINE fonte_neige_final
311!
312!****************************************************************************************
313!
314  SUBROUTINE fonte_neige_get_vars(pctsrf, fqcalving_out, &
315       fqfonte_out, ffonte_out)
316
317! Cumulate ffonte, fqfonte and fqcalving respectively for
318! all type of surfaces according to their fraction.
319!
320! This routine is called from physiq.F before histwrite.
321
322    INCLUDE "indicesol.h"
323!****************************************************************************************
324    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf
325
326    REAL, DIMENSION(klon), INTENT(OUT)      :: fqcalving_out
327    REAL, DIMENSION(klon), INTENT(OUT)      :: fqfonte_out
328    REAL, DIMENSION(klon), INTENT(OUT)      :: ffonte_out
329
330    INTEGER   :: nisurf
331!****************************************************************************************
332
333    ffonte_out(:)    = 0.0
334    fqfonte_out(:)   = 0.0
335    fqcalving_out(:) = 0.0
336
337    DO nisurf = 1, nbsrf
338       ffonte_out(:) = ffonte_out(:) + ffonte_global(:,nisurf)*pctsrf(:,nisurf)
339       fqfonte_out(:) = fqfonte_out(:) + fqfonte_global(:,nisurf)*pctsrf(:,nisurf)
340       fqcalving_out(:) = fqcalving_out(:) + fqcalving_global(:,nisurf)*pctsrf(:,nisurf)
341    ENDDO
342
343  END SUBROUTINE fonte_neige_get_vars
344!
345!****************************************************************************************
346!
347END MODULE fonte_neige_mod
348
349
350
Note: See TracBrowser for help on using the repository browser.